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ABSTRACT 


A FORTRAN IV computer program was written that gives the blade -to -blade solu- 
tion of the two-dimensional, subsonic, compressible {or incompressible), nonvlscous 
flow problem for a circular or straight infinite cascade of tandem or slotted turboma- 
chine blades. The blades may be fixed or rotating, The flow may be axial, radial, or 
mixed. The results include streamline coordinates, velocity magnitude and direction 
throughout the passage, and the blade -surface velocities. The method is based on the 
stream function using an iterative solution of nonlinear finite -difference equations. 
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FORTRAN PROGRAM FOR CALCULATING VELOCITIES AND STREAMLINES 
ON A BLADE-TO- BLADE STREAM SURFACE OF A 
TANDEM BLADE TURBOMACHINE 
by Theodore Katsanis and William D. McNally 
Lewis Research Center 

SUMMARY 

A FORTRAN IV computer program was written that gives the blade -to -blade solution 
of the two-dimensional, subsonic, compressible {or incompressible), nonviscous flow 
problem for a circular or straight infinite cascade of tandem or slotted turbomachine 
blades. The blades may be fixed or rotating. The flow may be axial, radial, or mixed, 
and there may be a change in stream -channel thickness in the through-flow direction. 

The program input consists of blade and stream -channel geometry, total flow condi- 
tions, inlet and outlet flow angles, blade -to -blade stream -channel weight flow, and the 
portion of this weight flow that passes between the front and rear tandem blades (through 
the slot). The output includes blade-surface velocities, velocity magnitude and direction 
at all interior mesh points in the blade -to -blade passage, and streamline coordinates 
throughout the passage. 

The method is based on the stream function. The simultaneous, nonlinear, finite - 
difference equations of the stream function are solved by using two major levels of iter- 
ation, The inner iteration consists of the solution of simultaneous linear equations by 
successive overrelaxation, using an estimated optimum over relaxation factor. The outer 
iteration then changes the coefficients of the simultaneous equations to correct for com- 
pressibility. 

This report includes the FORTRAN IV computer program with an explanation of the 
equations involved, the method of solution, and the calculation of velocities. Numerical 
examples are included to illustrate the use of the program, and to show the results which 
are obtained. 



INTRODUCTION 


An effort is being made to design compressors and turbines with smaller diameters, 
fewer stages, and fewer blades per stage. All these factors tend to increase diffusion. 
Therefore, it is desired to design blades with high diffusion, and at the same time to 
avoid flow separation. Several ideas for aerodynamic design to permit high diffusion 
without separation are being investigated, both theoretically and experimentally. Two 
promising concepts are the tandem blade and the slotted blade. 

In the design of tandem or slotted blade rows for compressors or turbines, an analy- 
sis is desirable which will give velocity distributions from blade to blade, and particu- 
larly over the blade surfaces, Stanitz (refs. 1 and 2) has shown that finite -difference so- 
lutions of the stream -function differential equation can be used to obtain these results. 
Computer programs have been written which generate coefficients for the difference equa- 
tions, solve the equations, and differentiate the resulting values of stream function to ob- 
tain velocities throughout the blade -to -blade passage and on the blade surfaces. This has 
been done previously by the first author (ref. 3) for a turbomachine with a single blade 
row without slots. 

This report extends the analysis of reference 3 to tandem or slotted blades. A com- 
puter program has been written to obtain the numerical solution for ideal, subsonic, com- 
pressible (or incompressible) flow for an axial-, radial-, or mixed-flow circular cascade 
of turbomachine blades. The program may also be used for a straight infinite cascade. 
The blades may be overlapping or nonoverlapping in the meridional flow direction and may 
be fixed or rotating. The program may also be used to analyze a turbomachine with one 
set of splitter blades (see section Mixed-Flow Impeller, p. 12). The coordinates used 
are meridional streamline distance and angular coordinate in radians. 

This report includes the FORTRAN IV computer program (called TANDEM) that was 
developed, with an explanation of the equations involved and the method of solution. A 
tandem axial gas-turbine rotor cascade and a mixed-flow impeller are analyzed to illus- 
trate the use of the program. The results obtained for the axial turbine are compared 
with experimental data. The impeller results are compared with previous analytical re- 
sults. The report is organized so that the engineer desiring to use the program needs to 
read only the sections MATHEMATICAL ANALYSIS, NUMERICAL EXAMPLES, and 
DESCRIPTION OF INPUT AND OUTPUT. Information of interest to a programmer is 
contained in the sections DESCRIPTION OF INPUT AND OUTPUT and PROGRAM PRO- 
CEDURE and in the appendixes. 

A TANDEM source deck on tape is available from COSMIC (Computer Software 
Management and Information Center), Computer Center, University of Georgia, 

Athens, Georgia 30601. The program number can be obtained from the authors. 
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SYMBOLS 


A 


a ij 


a Q> a l‘ a 2’ a 3’ 
a 4 J a 12 > a 34 


} 


h 

b 12’ b 34 

C P 

h 


k 


m 

n 

R 

r 


s 

T 

u 


u m 


coefficient matrix, eq, (A7) 
typical element of matrix A 

coefficients in eq. (A2) 


stream-channel thickness normal to meridional streamline, meters 
quantities in eq. (A2) 

specific heat at constant pressure, joule/(kg)(°K) 

spacing between adjacent points, eqs. (Al) to (A4); see fig. 17 



constant vector, 


, eq. (A7) 



meridional streamline distance, see figs. 2 and 3 
number of unknown mesh points 
gas constant, joule/{kg)(°K) 

radius from axis of rotation to meridional stream -channel mean line, 
meters 

angular blade spacing or pitch, rad 
temperature, °K 
stream function 



discrete approximation to stream function at n mesh points, 


. t T ) 
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V absolute fluid velocity, meters/sec 

W fluid velocity relative to blade, meters/sec 

w mass flow per blade flowing through stream channel, kg/sec 

z axial coordinate, meters 

a angle between meridional streamline and axis of rotation, rad; see fig. 1 

0 angle between relative velocity vector and meridional plane, rad; see fig. 1 

y specific -heat ratio 

n outer normal to region 

6 relative angular coordinate, rad; see fig 1 

9 

X prerotation (rV,,) , meters /sec 

p'in 

p density, kg/meters 

£2 over relaxation factor, eq. (A8) 

u) rotational speed, rad/sec; see fig. 1 

Sub scripts: 

cr critical velocity 

1 dummy variable 

in inlet or upstream 

j dummy variable 

le leading edge 

m component in direction of meridional streamline 

out outlet or downstream 

te trailing edge 

6 tangential component 

0, 1,2, 3, 4 quantities at these locations in finite -difference expression, eqs. (Al) to 
(A6); see fig. 17 

Superscripts: 

T transpose of vector or matrix 

' absolute stagnation condition 

" relative stagnation condition 
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MATHEMATICAL ANALYSIS 

It is desired to determine the flow distribution through a stationary or rotating cas- 
cade of tandem blades on a blade -to-blade surface. The following simplifying assump- 
tions are used in deriving the equations and in obtaining a solution: 

(1) The flow is steady relative to the blade. 

(2) The fluid is a perfect gas or is incompressible. 

(3) The fluid is nonviscous. 

(4) There is no loss of energy. 

(5) The flow is absolutely ir rotational. 

(6) The blade -to -blade surface is a surface of revolution. (This does not exclude 
straight infinite cascades. ) 

(7) The velocity component normal to the blade -to -blade surface is zero. 

(8) The stagnation temperature is uniform across the inlet. 

(9) The velocity magnitude and direction is uniform across both the upstream and 
downstream boundaries. 

(10) The relative velocity is subsonic everywhere. 
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The flow may be axial, radial, ♦ or mixed, and there may be a variation in the stream- 
channel thickness b in the through -flow direction. The proportion of flow between the 
front and rear blades must be specified as an input to the program. This input may be 
difficult for the user to estimate; however, correlation with experimental work may yield 
more reliable values. 

The coordinate system is shown in figure 1. Since the variables r and z are not 
independent on a stream surface, one variable can be eliminated. Therefore, r and 0 
or z and 0 could be used as independent variables. However, for generality, it is 
better to use the meridional streamline distance m in place of r and z as an indepen- 
dent variable (see fig. 2). Then, m and 9 are the two basic independent variables. A 
stream channel is therefore defined by specifying a meridional streamline radius r and 
a stream -channel thickness b at several meridional locations m (see fig. 3). 

For the mathematical formulation of the problem, the stream function is used. The 
stream function u used herein is related to the stream function ^ defined in refer- 
ence 4 by u = -i/'/w. With this substitution in equation 12(9) of reference 4 we obtained 
the basic differential equation which must be satisfied by the stream function under the 
given assumptions: 


0ja + 3 2 u 
r 2 30 2 3m 2 


1 dp 3u 
r 2 p 30 30 


sin a _ J_ 3(bp) 
r bp 3m 


3u . 2bpq> 
3m w 


sin a 


( 1 ) 


The stream function u has the value 0 on the upper surface of the leading blade and 1 on 
the lower surface of the leading blade. Also, the derivatives of the stream function sat- 
isfy the equations 


w fl 

3m w w 

(2) 

3u _ bpr y, 

77 m 

30 w 

(3) 


For the solution of equation (1), a finite region is considered (as indicated in fig. 4) 
with the condition that the flow along corresponding upper and lower portions of the bound- 
ary is the same. For example, the flow along AB is the same as along NM* Also, it is 
assumed that AN is sufficiently far upstream so that the flow is uniform along this bound- 
ary, and that the flow angle is known. Similarly, it is assumed that the flow is uni- 
form along GH, and that the flow angle /3 Qut is known. For an actual blade row, /3 ^ 

may usually be determined by means of experimentally determined rules. Also, it is as- 
sumed the flow split is known; that is, the percentage of flow which passes between the 
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lb* Nonoverfap case. 
Figured * Typical finite flow region. 


front and rear blades. Specifying /3 t and the flow split is mathematically equivalent 
to specifying the locations of the stagnation points on the trailing edges of both blades. 

Since equation (1) is elliptic for subsonic flow, boundary conditions for the entire 
boundary ABCDEFGHUKLMNA are required. Along BC, u = 0; along LM, u = 1; along 
EF, u is equal to the negative of the fraction of weight flow through JKL; and along JJ, 
u is equal to the fraction of weight flow crossing a line joining C and J. Along AB, CD, 
FG, HI, KL, and MN, a periodic condition exists; that is, the value of u along MN, KL, 
and HI is exactly 1. 0 greater than it is along AB, CD, and FG. The same condition holds 
along DE and JK in the nonoverlapping case (fig. 4(b)). 

Along AN and GH, 3u/ dr) is known, where rj is in the direction of the outer normal. 
From equations (2) and (3), since W^/ W m = tan /3, 

JH. = . tan (3 
3m r 36 


I 
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Along AN and GH, 


where s is the angular blade 


3u _ u(N) - u(A) _ l 
36 s s 

spacing, so that 



along AN 


along GH 


( 4 ) 

(5) 


These are the boundary conditions required to determine a solution to equation (1). The 
method used for the numerical solution of equation (1) is described in appendix A. The 
numerical solution involves two levels of iteration because equation (1) is nonlinear. The 
inner iteration is required to solve equation (1) when it is linearized, and the nonlinear 
solution is approached by the outer iteration. 

After computing a numerical solution to equation (1) in a given flow region, the ve- 
locity at any point can be computed from equations (2) and (3) by using numerical differ- 
entiation. The streamlines are located by the contours of equal stream -function values. 


NUMERICAL EXAMPLES 

To illustrate the use of the program and the type of results which can be obtained, 
two numerical examples are given. The first example is an axial-flow turbine and the 
other is a mixed -flow impeller. 


Axial- Flow Turbine Rotor Cascade 

This example is a two-dimensional axial-flow turbine cascade currently undergoing 
testing at Lewis Research Center. This blade is a modified version of a tandem blade 
reported in reference 5. It has a blunt leading edge on the rear blade in order to achieve 
a converging channel between the blades, and it has a wider slot than that reported in 
reference 5. 

The blade shape in m, e coordinates and the blade -to -blade solution region are 
shown in figure 5. Input for this example is given in table I. Blade -surface velocities 
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TABLE L - INPUT FOR AXIAL- FLOW TURBINE ROTOR CASCADE 


MODIFIED TANDEM AXIAL TURBINE ROTOR 
SAM Aft TIP 

i.4000000 207,05300 238*15000 

BETA 1 8 E TAG CHOftQF 

48.000000 -A 7,000000 0 . 28470006-0 t 

MB I MBO MB 1 2 MB02 MM MBS I NftL NRSP 
10 32 29 49 58 20 7f» 2 


RHOIP 

1.2250000 

STGRF 

0.213 33006-01 


WTFL 

0, 31520006-01 
CHGRDft 

0.25150006-01 


WTFLSP 

0* i 1347006-0 l 
STGRR 

-0.54590006-01 


QM6GA 

-0 

ML6R 

0.24410006-01 


QRf 

0 

THLER 

-O.3&07Q006-02 


BLADE SURFACE 1 — UPPER SURFACE - FRONT BLADE 

ftil ftOt 8ETU B6T01 

0,76200006-03 0. 38 100006-03 50.000000 -29.400000 

MSP1 ARRAY 

-0 0, 25700006-02 0*76500006-02 0. 1527000E-Q 1 

THSP1 ARRAY 


-0 


0,92500006-0 2 0. 21180006-01 


0*29380006-01 


SPLN01 

7.0000000 

0,20350006-01 0,25430006-01 -0 

0.30200006-01 0.26430006-01 -0 


BLAD6 SURFACE 2 — LOWER SURFACE - FRONT BLADE 

RI2 R02 86 T 12 BE102 

0*76200006*03 0,38100006-03 25.000000 -6.9000000 

MSP2 ARRAY 

*0 0.76500006-02 0,20350006-01 0,25430006-01 

THSP2 ARRAY 


-0 


0*71400006-02 0,20390006-01 


0. 20940006-01 


SPLN02 

5,0000000 

-0 

-0 


BLADE SURFACE 3 — UPPER SURFACE - REAR BLADE 


R I 3 

403 

BET 13 

B6T03 

0. 1 7780006*02 

0.3810000E-03 

-8*1000000 

-48.800000 

MSP3 ARRAY 

0 

0.61000006-02 

0. 16260006-01 

0 

THSP3 ARRAY 

0 

0.16450006-02 

-0,24630006-01 

0 


SPLN03 

4.0000000 


BLADE SURFACE 4 — LOWER SURFACE - REAR BLADE 


414 

R04 

66TI4 

86T04 

0.17780006-02 

0.38100006-03 

-19.700000 

-42.500000 

HSP4 ARRAY 

0 

0,61000006-02 

0.13720006-01 

0 

THSP4 ARRAY 

0 

-0.12000006-01 

-0.27450006-01 

O 


SPLN04 

4,0000000 


MR ARRAY 

■ 1,0000000 

RHSP ARRAY 
0*3238500 

BESP ARRAY 
0*10000006-01 


1.0000000 
0, 3238500 
0.10000006-01 


BLOAT AAJ40K ERSOR STftPN SLCftD INTVL SURVL 

1 1 2 2 2 2 3 



Ill 



figure 6, - Surface velocities on tandem axial turbine blade. 


are plotted in figure 6, where comparison is made with unreported experimental data for 
the Lewis turbine cascade* There is close agreement between computed and experimental 
values on all four blade surfaces. 

Execution time was 10 minutes for this example, and it required 16 outer iterations 
for final convergence to the compressible solution. 


Mixed-Flow Impeller 

This example is taken from reference 6. In reference 6 a similar stream -function 
analysis was made. The mesh was set up graphically, and the coefficients were calcu- 
lated by hand. The solution of the finite-difference equations was obtained by relaxation 
on a computer. The analysis was done on a blade -to -blade surface of revolution midway 
between hub and shroud. 

The coordinates of the stream channel and the stream -channel radial thickness are 
given by equations (!) and (2) of reference 7. The radial stream -channel thickness was 
corrected to obtain the normal thickness required by this program. The hub-shroud pro- 
file is shown in figure 7. The blade shape and mesh arrangement are shown in figure 8. 
Input for this example is given in table n. In figure 9 the blade-surface velocities ob- 
tained by the TANDEM program are compared with those obtained originally in refer - 
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Figure 7. - Hub-shroud profile of mixed-flow impeller showing meridional section of steam tube. 
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TABLE II. - INPUT FOE MIXED-FLOW IMPELLER 


MIXED FLOW 

impeller inasa 

TN U-U&6* 






GAM 

AR 

TIP 

RHOIP 

WIFL 

WTFLSP 

OMEGA 

QRF 

1*5000000 

1000*0000 

1000000.0 

1,0000000 

0.30420OO6-02 

0.13516006-02 

796*00000 

0 

BETA! 

8EIA0 

CHORD? 

STGRF 

CHOROR 

STGRR 

ML6R 

THLER 

-B4.8BQOOO 

MB I M0O mi? 
10 47 28 

-43.000000 
MBQ2 MM MSB 1 
47 57 28 

0.1055500 
McSL NR5P 
8 18 

-2*6290000 

0,5664000 E-Q l 

-0,6649000 

0, 4891000E-01 

-2.3434000 


blade SURFACE 

1 — UPPER SURFACE - FRONT 

BLADE 



Rll 

RGt 

BETH 

BET01 

SPLN01 


0 * 9 140000E-01 

0*18460006-02 

-80.000000 

-49*000000 

6.0000000 


MSP 1 ARRAY 

0 

0*12140006-01 

0.2651000E- 

01 0,47660006-01 

0.73600006-01 

0 

THSPl ARRAY 

0 

-O.625O0QD 

- 1 ,2330000 

- i * 8182000 

-2,2750000 

0 


GLADE SURFACE 

2 LOWER SURFACE - FRONT 

8LA0F 



ft I Z 

R02 

BET 1 2 

86 T02 

SPLN02 


0*91400006-03 

0*18460006-02 

-83 * 000000 

-41.500000 

6*0000000 


MSP 2 ARRAY 






0 

0-78800006-02 

0* 2004000E- 

01 0*40060006-01 

0*68280006-01 

0 

THSP2 ARRAY 






0 

-0*6310000 

- 1 * 23 10000 

-1.8206000 

-2*2954000 

0 


blade Surface ^ — upper surface - rear 

blade 



R I 3 R03 

BET 1 8 

B6T03 

SPLND3 


0. 1 3280006-02 0* 1 7530006-02 

-60.500000 

-51,500000 

6*0000000 


MSP3 ARRAY 

0 0,130 70OQE-Q1 

0.2552G0GE 

-01 0,41720006-01 

0 * 5280000E-0 1 

0 

THSP3 ARRAY 

0 -0*t670000 

-0*3370000 

-0,5262000 

-0 * 6269000 

0 


BLADE SURFACE 4 -- LOWER SURFACE - REAR *L 

RI4 RQ4 B6TI4 

0*1 32SOOOE -02 0*1 7530Q0E-02 -60*000000 

M$P4 ARAAV 

0 0* 107 30006*0 1 0*24930006-01 


ADE 

BE TO A 

- 40*500000 


O,41720QQE«O1 


THSP4 ARRAY 


-0,2010000 

-0*4070000 


-0*5315000 


SPLNQ4 

5*0000000 

0 

0 


MR ARRAY 
o*3i240ooe-oi 
0,51150006-01 
0,1272600 

RM5P ARRAY 
0, 75060006-0 I 
0* 944 7000E-0 i 
0* L 360200 

BE S P ARRAY 
0.|0533OOt“0t 
0 * 32350OQE -02 
0.825OQOOE -03 


-0, 15 1 4000E-0 l 
0,59640006-01 
0,1407300 

Q.76620OQE-01 

0.98200006-01 

0*1445700 

0,10045006-01 
0*2 7280Q0E -02 
0.72400OOE-G3 


0 . 2500000E -03 
0 * 682B000E -01 


o*7S74onoe-oi 

0,1022800 


0 • 8 72 4000 E -02 
Q*2299QO0E-02 


O.106500QE-O1 

O.77Q900OEKU 


0 , 80910006-01 
0.1067400 


0-7420000E-02 

0,19360006-02 


0*18530006-01 

O-86Q70OO6-01 


0*62940006-01 

0*1114600 


0,63160006-02 

G.1629Q0QE-02 


0.265 1000E- 01 
0.9524000E-01 


0,85310006-01 

0.1165600 


O.53540OQE-O2 

0*13700006-02 


0,34600Q0&~0i 

0.1046100 


0,88020006-01 

0, 1220000 


0*45320006-02 
0*1 15 I0006-G2 


o. 42 aioooe~oi 

0*1141700 


0.91080001-01 
0. I 277800 


0*38310006-02 

0.97900006-03 


Vi 


bloat aanok ERSrw STRFn SLCRO imtvl 

l 1 2 2 2 2 


SURVL 

3 


TANDEM program 
Ref. 6 


Slade surface 2 



I I I I I 

.4 .6 

Fraction of distance between stagnation points 

(a) Full blade. 

Figure 9. - Velocity distribution for mixed -flow pump impeller. 








ence 6. There is good agreement over most of the blade. Minor discrepancies are prob- 
ably due to slight differences in the boundary conditions (weight flow split and downstream 
flow angle). The heights of the peaks near the leading edges are uncertain because the 
radii are small compared to the mesh spacing. 

Execution time was 2 minutes for this example. It required only one outer iteration, 
since flow was incompressible. 


DESCRIPTION OF INPUT AND OUTPUT 

The computer program requires as input a geometrical description in m,Q coordi- 
nates of the tandem blade segments, a description in m, r coordinates of the stream 
channel through the blades, appropriate gas constants, and operating conditions such as 
inlet temperature and density, inlet and outlet flow angles, weight flow, and rotational 
speed. An estimate of the portion of the weight flow which passes between the tandem 
blades must also be given. Output obtained from the program includes velocity magnitude 
and direction at all interior mesh points in the blade -to-blade passage, blade -surface 
velocities, stream -function values throughout the blade-to-blade region of solution, and 
streamline locations. 


Input 


Figure 10 shows the input variables as they are punched on the data cards. The 
first input data card is for a title, which will serve for problem identification. The 
remaining cards are for input variables. There are two types of variables, geometric 
and nongeometric. The geometric input variables are shown in figures 11 to 13. Fur- 
ther explanation of the input variables is given in the section Instructions for Prepar- 
ing Input. 

The input variables are as follows: 


GAM specific -heat ratio, y 

AR gas constant, joule/(kg)(°K) 

TIP inlet total temperature, T^, °K 

RHOIP inlet total density, p! Q , kg/meter 

WTFL mass flow per blade for stream channel, kg/sec 

WTFLSP portion of stream -channel mass flow per blade which flows across 

the boundary JKL between the front and rear blades, kg/sec; see 
fig. 11 
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60j61 70171 Ml 


WFLSP 

STCRR 





Figure 10. - Input form. 
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CHORD? 


STGRRH 


U — CHORDR *1 L 

- m 


(al Overlapping case. 



Ibi Nonovertapping case. 

- Geometric input variables for blade -to -blade flow region. 





Figure 12, - Geometric input variables on blade. Angles BETH, 2, 3, 4 and 8ET01, 2, 3, 4 
must be given as true angle JJ, not as angles measured in m, 9 plane. Either use 
tang' r d9/dm to obtain fi, or measure the true angle. 



OMEGA 

ORF 

BETAI 

BETAO 

CHORDF 


rotational speed, co t rad/ sec (Note that w is negative if rotation is 
in the opposite direction of that shown in fig. 1. ) 

value of overrelaxation factor S2 to be used in equation (A8) (If 
ORF = 0, the program calculates an estimated value for the over- 
relaxation factor; see p. 25 and appendix A for discussion. ) 

inlet flow angle 0, along BM with respect to m -direction, deg; see 
fig. 11 

outlet flow angle *te along FI with respect to m -direction, deg; see 
fig. 11 

overall length of front blade in m -direction, meters; see fig. 11 
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STGRF 


CHORDR 

STGRR 


MLER 

THLER 

MBI 

MBO 

MBI2 

MB02 

MM 

NBBI 

NBL 

NRSF 

RI1, RI2, 

RI3, RI4 

ROl, R02, 
R03, R04 

BETH, BETI2, 
BETI3,BETI4 


BETOl, BET02, 
BET03, BET04 


angular 6 -coordinate for center of trailing-edge circle of front blade 
with respect to center of leading-edge circle of front blade, rad; 
see fig. 11 

overall length of rear blade in m -direction, meters; see fig. 11 

angular 9 -coordinate for center of trailing-edge circle of rear blade 
with respect to center of leading-edge circle of rear blade, rad; 
see fig. 11 

m -coordinate of leading edge of rear blade with respect to leading 
edge of front blade, meters; see fig. 11 

angular 6 -coordinate of leading edge of rear blade with respect to 
leading edge of front blade, rad; see fig, 11 

number of vertical mesh lines from AN to BM inclusive; see fig. 11 

number of vertical mesh lines from AN to CL inclusive; see fig. 11 

number of vertical mesh lines from AN to EJ inclusive; see fig. 11 

number of vertical mesh lines from AN to FI inclusive; see fig. 11 

total number of vertical mesh lines in m -direction from AN to GH, 
maximum of 100; see fig. 11 

number of mesh spaces in 9 -direction between AB and MN, maxi- 
mum of 50; see fig. 11 

number of blades 

number of spline points for stream -channel radius (RMSP) and thick- 
ness (BESP) coordinates, maximum of 50; see fig. 13 

leading-edge radii of the four blade surfaces, meters; see fig, 12 

trailing-edge radii of the four blade surfaces, meters; see fig. 12 

angles (with respect to m -direction) at tangent points of leading-edge 
radii with the four blade surfaces, deg; see fig. 12 (These must 
be true angles in degrees. If angles (i. e. , de/dm) are measured 
in the m,e plane, BETI1,2, 3, 4 can be obtained from the relation 
tan /3 = r 69 / dm. ) 

angles (with respect to m -direction) at tangent points of trailing-edge 
radii with the four blade surfaces, deg; see fig. 12 (These must 
also be true angles in degrees, like BETX1,2, 3, 4.) 
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SPLNOl, SPLN02, number of blade spline points given for each surface as input, maxi- 
SPLN03, SPLN04 mum of 50 (These include the first and last points (dummies) that 

are tangent to the leading- and trailing-edge radii (fig. 12). ) 


MSP1, MSP2, 
MSP3, MSP4 


THSP1, THSP2, 
THSP3, THSP4 


arrays of m -coordinates of spline points on the four blade surfaces, 
measured from blade leading edges, meters; see fig. 12 (The 
first and last points in each of these arrays can be blank or have a 
dummy value, since these points are calculated by the program. 

If blanks are used, and the last point is on a new card, a blank 
card must be used. ) 

arrays of 6 -coordinates of spline points corresponding to MSP1, 
MSP2, etc., rad; see fig. 12 (Dummy values are also used in 
positions corresponding to those in MSP1, MSP2, etc. ) 


MR 


RMSP 

BESP 


array of m -coordinates of spline points for stream -channel radii 
and stream -channel thicknesses, meters; see fig. 13 (MR is 
measured from leading edge of front blade. These coordinates 
should cover the entire distance from AN to GH and may extend 
beyond these bounds. The total number of points is NRSP. ) 

array of r -coordinates of spline points for stream -channel mean 
streamline, corresponding to the MR array, meters; see fig. 13 

array of stream -channel normal thicknesses corresponding to the 
MR and RMSP arrays, meters; see fig. 13 


The remaining variables, starting with BLDAT, are used to indicate what output is 
desired. A value of 0 for any of these variables will cause the output associated with that 
variable to be omitted. A value of 1 will cause the corresponding output to be printed for 
the final outer iteration only; a value of 2, for the first and final iterations; and a value 
of 3, for all outer iterations. Care should be used not to call for more output than is 
really useful. The following list gives the output associated with each of these variables: 

BLDAT all geometrical information which does not change from iteration to iteration; 
i. e. , coordinates and first and second derivatives of all blade- surface 
spline points; blade coordinates and blade slopes where vertical mesh lines 
meet each blade surface; radii and stream -channel thicknesses corresponding 
to each vertical mesh line; m -coordinate, stream -channel radius and thick- 
ness, blade- surface angles and slopes where horizontal mesh lines intersect 
each blade; and ITV and IV arrays (internal variables describing the location 
of the blade surfaces with respect to the finite-difference grid) 
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AANDK coefficient array A, vector k, and indexes of all adjacent points for each 

point in the finite -difference mesh (This information is needed only for de- 
bugging the program. ) 

ERSOR maximum change in stream function at any point for each iteration of SOR 
equation (eq. (A8)) 

STRFN value of stream function at each unknown mesh point in region 

SLCRD streamline 9 -coordinates at each vertical mesh line, and streamline plot 

INTVL velocity and flow angle at each interior mesh point 

SURVL m -coordinate, surface velocity, flow angle, distance along surface, and 

W/W based on meridional velocity components where each vertical mesh 

V* 

line meets each blade surface; m-coordinate, surface velocity, flow angle, 
distance along surface, and W/W cr based on tangential velocity components 
where each horizontal mesh line meets each blade surface; plot of blade- 
surface velocities against meridional streamline distance, m (It is suggested 
that SURVL = 3 be used. This will give surface velocities after each outer 
iteration, so that satisfactory velocities may be obtained even when final con- 
vergence is not reached. ) 

Instructions for Preparing input 

Units of m easurement. - The International System of Units (ref. 8) is used through- 
out this report. However, the program does not use any constants which depend on the 
system of units being used. Therefore, any consistent set of units may be used in pre- 
paring input for the program. For example, if force, length, temperature, and time are 
chosen independently, mass units are obtained from force equals mass times accelera- 
tion. The gas constant R must then have the units of force times length divided by mass 
times temperature (energy per unit mass per deg). Density is mass per unit volume, and 
weight flow is mass per unit time. Output then gives velocity in the chosen units of length 
per unit time. Since any consistent set of units can be employed, the output is not labeled 
with any units. 

Blade and stream -ch annel geometry. - The upper and lower surfaces of the front and 
rear tandem blades are each defined by specifying three things: leading- and trailing-edge 
radii, angles at which these radii are tangent to the blade surfaces, and m- and 6 - 
coordinates of several points along each surface. These angles and coordinates are used 
to define a cubic spline curve fit (ref. 9) to the surface. The standard sign convention is 
used for angles, as indicated in figure 12. 
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A cubic spline curve is a piecewise cubic polynomial which expresses mathematically 
the shape taken by an idealized spline passing through the given points. Reference 9 de- 
scribes a method for determining the equation of the spline curve. When this method is 
used, few points are required to specify most blade shapes accurately, usually no more 
than five or six, in addition to the two end points. As a guide, enough points should be 
specified so that a physical spline passing through these points would accurately follow 
the blade shape. This means that the spline points should be closer where there is large 
curvature and farther apart where there is small curvature. 

The coordinates for either surface of a particular blade segment are given with re- 
spect to the leading edge of that segment, the leading edge being defined as the furthest 
point upstream on the blade segment. 

The mean stream surface of revolution (as seen in the meridional plane, fig. 13) and 
the stream -channel thickness are also fitted with cubic spline curves. The m -coordinates 
for the mean stream surface are independent of the m -coordinates for blade surfaces. 

Inlet and outlet flow angles . - The values of j3 le and (3^ e are given as average val- 
ues on BM and FI, respectively. If the flow is axial, these flow angles are the same as 
the flow angles at AN and GH. If flow is radial or mixed, and these angles are not known 
on BM and FI, /3^ e and /3^ g must be calculated by equation (B14). 

Defining mesh . - A finite -difference mesh is used for the solution of equation (1). A 
typical mesh pattern (that used in example 1) is shown in figure 14. The mesh spacing 
and the extent of the upstream and downstream regions are determined by the values of 
MBI, MBO, MBI2, MB02, and MM of the input (fig- 10)- The mesh spacing must be cho- 
sen so that there are not more than 2000 unknown mesh points. 

Values of MBI, MBO 2, etc. , should be determined so that the mesh which results has 
blocks which are approximately square. To achieve this, a value for NBBI is first chosen 
arbitrarily (15 to 20 is typical). NBBI is the number of mesh spaces spanning the blade 
pitch s, where s = 27r/NBL. Dividing s by NBBI gives the mesh spacing HT in the 8- 
direction in radians. Multiplying HT by an average radius (RMSP) of the stream channel 
gives an average value for the actual mesh spacing in the Q -direction. CHORDF, 
CHORDR, and MLER should then be used with this tangential mesh spacing to calculate 
the approximate number of mesh spaces in the various regions along the meridional axis. 
This will give MBO, MB 12, and MB02, once MBI is chosen. Generally, MBI is given a 
value of 10; MM, likewise, is usually given a value 10 more than MBO 2. 

Over re taxa tion factor . - ORF is the relaxation factor used in each inner iteration in 
the solution of the simultaneous finite -difference equations (A7). ORF may be set equal 
to 0, or to some value between 1 and 2. ORF is usually given as 0 for the initial run of a 
given blade geometry and mesh spacing (MBI, NBBI, etc.). In this case the program 
uses extra time and calculates an optimum value for ORF. It does this by means of an 
iterative process, and on each iteration the current estimate of the optimum value for 
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ORF is printed. The final estimate is the one used by the program for ORF. If the user 
does not change the mesh indexes MB I, MBO, MB 12, MB02, MM, and NBBI between runs, 
even though blade geometry or other input does change, he may use this final estimate of 
ORF in the input, saving the time used in its computation. In all cases, if ORF is not 0, 
it should have a value greater than .1 and less than 2. 

Actually, the value of ORF is not as critical as the user might think. It gets more 
critical as the optimum value gets close to 2. For any run of a given set of data, only 
small changes will occur in the rate of convergence in SOR as long as the d if ference 
2.0 - ORF is within 10 percent of its optimum value. A further theoretical discussion 
of the overrelaxation factor is presented in reference 11 (p. 78). 

Format for input data . - All the numbers on the card beginning with MBI and on the 
card beginning with BLDAT are integers (no decimal point) in a five-column field {see 
fig. 10). These must all be right adjusted. The input variables on all other data cards 
are real numbers (punch decimal point) in a ten -column field. 

Incompressible flow . - While the program is written for compressible flow, it can 
be easily used for incompressible flow. To do so, specify GAM =1.5, AR = 1000, and 
TIP = 10° as input o This results in a single outer iteration of the program to obtain the 
stream-function solution. 

Straight infinite cascade , - The program is as easily applied to straight infinite cas- 
cades as to circular cascades* Since the radius and number of blades (NBL) for such a 
cascade would actually be infinite, an artificial convention must be adopted. The user 
should pick a value for NBL, for instance 20 or 30* Then, since the blade pitch equal to 
sr is known, an artificial radius can be computed from 

r = NBL*(sr) 

2n 

This r should be used to compute the 6 -coordinates required as input (THSP1, . . . , 
THSP4, STGR1, 3TGR2, THLE2) by dividing coordinates in the tangential direction by 

r* 

Axial flow . - For a two-dimensional cascade with constant stream -channel thickness, 
constant values should be given for the RMSP and BESP arrays* Only two points 
are required for each of these arrays in this case* The two values of MR should be 
chosen so that they are further upstream and downstream than the boundaries AN and GH* 
The two values of RMSP and BESP should equal the constants r and b. 


Output 


Sample output is given in table III for the axial -flow turbine example. Since the com- 
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plete output would be lengthy, only the first few lines of each section of output are repro- 
duced herein. Most of the output is optional and is controlled by the final input card, as 
already described. In many instances output labels are simply internal variable names 
which are defined in the Main Dictionary. 

Each section of the sample output in table HE has been numbered to correspond to the 
following description: 

(1) The first output is a listing of the input data. All items are labeled as on the in- 
put form (fig. 10). 

(2) This is the output corresponding to BLDAT* (See the list of input variables 
and the Main Dictionary for variable name definitions. ) 

(3) The relative free -stream velocity W, the relative critical velocity W cr , and 

the maximum value of the mass flow parameter pW (corresponding to W = W cr ) are 
given at the leading edge of the front blade (BM) and the trailing edge of the rear blade 
(FI). The inlet (outlet) free -stream flow angle P- n (|3 j.) at boundary AN (GH) is given. 

These angles are based on the input angles BETAI (0 le ) and BETAO 03 te ). 

(4) These are calculated program constants, including the pitch from blade to blade, 
the mesh spacing in all solution regions, the minimum and maximum values of IT in the 
solution region (ITMIN and ITMAX), and the value of the prerotation A (eq. (B8)). 

(5) This is the number of mesh points in the entire solution region at which the 
stream function is unknown. 

(6) This is the boundary value (BV) of the stream function on each of the four blade 
surfaces. 

(7) This is the output corresponding to AANDK. 

(8) If the program calculates an optimum over relaxation factor 12 (i. e. , ORF = 0 in 
the input), the successive estimates to the optimum value of ORF are printed. The last 
printed value of the estimated optimum ORF is the value of 12 (ORF) used by the program. 

(9) This is the output corresponding to ERSOR. 

(10) This is the output corresponding to STRFN. 

(11) This is the total execution time after obtaining the stream -function solution 
for each outer iteration. 

(12) This is the output corresponding to SLCRD. 

(13) This is the output corresponding to INTVL. 

(14) This gives the maximum relative change in the density p for each outer itera- 
tion. 

(15) This is the output corresponding to SURVL. 

(16) This is the total execution time after all calculations are completed for an outer 
iteration. 
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TABLE m. - SAMPLE OUTPUT 


f EDIFIED TANDEM AXIAL TURBINE 

ROTOR 






SAM 

AR 

TIP 

RHOIP 

WTFL 

WTFLSP 

OMEGA 

ORF 

1,4000000 

287,05300 

268,15000 

1,2250000 

0, 3152000E-0 1 

0, 1 I34700E-01 

-0 

0 

SETAI 

REIAT 

chordf 

STGRF 

CHORDR 

STGRR 

MLER 

THLER 

40,000000 

-47,000000 

0,2fl47000E-01 

Q-2 L33300E-01 

0,25150006-01 

-0,54590006-01 

0,244 L000E-01 

-0.360700DE-02 

MBI MBO M812 MB02 MM N08 [ NBL 

NRSP 






10 32 29 49 58 20 76 

2 






BLADE SURE ACE 

1 — UPPER SURFACE - FRONT BLADE 





R I 1 

RD1 

BETI L 

BETD1 

SPtNOL 




0.76200006-03 

0- 381GO0OE-O3 

50,000000 

-2 9* 400000 

7,0000000 




MSP1 ARRAY 








-o 

0-2570000E-G2 

0.765QQO0E-02 

0, 1527000E-01 

0-2O35OOOE-OI 

0,2S43000E-0k 

-0 


IH5P1 ARRAY 








-0 

G.9250000E-02 

0,21 laoooc-ot 

0,29860006—0 1 

0, 302000GE-01 

0, 2643GOOE-Q1 

-0 


BLADE SURFACE 

2 -- LOWER SURFACE - FRONT BLADE 





RI2 

RQ2 

BET12 

6ET02 

SPLN02 




0, 7620000E-Q3 

0, 3B10000E-Q3 

25,000000 

-6,9000000 

5-0000000 




MSP 2 ARRAY 








-0 

0, 765000QE-02 

0.2Q35000E-01 

0,25430006-01 

-0 




THSPZ ARRAY 








! -0 

0* 7 14000QE-02 

0,20390006-01 

0-2G94000E-01 

-0 





BLADE SURFACE 

3 — UPPER SURFACE - REAR SLADE 


RI3 

R03 

BETI3 8ETD3 

SPLN03 

0, L778000C-02 

0,30100006-03 

-0,1000000 -4B, 800000 

4.0000000 

HSP3 ARRAY 




0 

0, 6 LOOOOOE-02 

0,16260006-01 0 


IHSP3 ARRAY 




0 

0-1640000E-02 

-0,24630006-01 0 



glade surface 

4 — LJWER SURFACE - REAR RLAOE 


RI4 

R04 

SET I 4 

SETQ4 

SPLN04 

0.1778000E-02 

O.301OOOOE-O3 

-L9, 700000 

42,500000 

4,0000000 

MSP4 ARRAY 





0 

0,61000006-02 

0. 1372000E-01 

0 


THSP4 ARRAY 





0 

-0, 1200000E-01 

-0,27450006-01 

0 



MR ARRAY 

-1,0000000 1* ooooooo 

RMSP ARRAY 

0,3230500 0,323fl500 

BESP ARRAY 

o,looooooE-oi 0 - looooooe-oi 



BLOAT 

AANDK 

EftSOR 

STRF N 

SLCRD 

1NTVL 

SURVL 

l 

1 

1 

2 

l 

2 

2 

3 


TABLE IE. -Continued. SAMPLE OUTPUT 


BLADE DATA, AT INPUT SPLICE POINTS 


O*17027F-O3 
0,257QOE-02 
0-765OOE-O2 
0* 1 52 7 OP- 01 
0-Z0350E-0L 
0*254 30F-01 
0*2B276E-0l 


BLADE SURF 
THETA 

0, 1 5 l 24E -02 
O,9250OE“Q2 
0*21 190E -0 L 
0. 29800E-D i 
0* 30200E-0 1 
0.26430F-01 
0.22358E-QI 


C£ l 
DERIVATIVE 
3*67996 
2-00 161 
L- 0 3901 
0-47565 
“0.33906 
- 1*15630 
-1*73991 


2ND DER IV. 
-448* 310 
-219*276 
-191* L99 
-166*633 
- 159.115 
- 167*823 
-241*946 


M 

0*lOfl4OF“O2 
0-76500E-02 
2 v 0.20350E-01 
] O.Z5430F-01 

0* 28Q43E-01 


0LADE SURFACE 2 
THETA DERIVATIVE 

0.80541E-01 1,43939 

O*»98I3E-01 1.38132 

0.10306 0,43322 

0,10361 “0. I 8877 

0. 10289 -0.37367 


2 MO DERI V 

-7.45650 

“10-3836 

“138*924 

-105-953 

-35*5599 


M 

0.26439E-01 

0*305108-01 

0.4O67QF-O1 

0*49466F“Gl 


BLADE 
THETA 

0* 1B284E-02 
“0- I9670E-02 
-0-28237E-01 
*0.574228-01 


SURFACE 3 

derivative 

-0*43947 
- 1-49683 
“3. 17480 
-3.52722 


2ND DERI V * 

-206*721 

-3L2-6B1 

-17,6280 

-62,5065 


M 

0.25589E-01 
0-305 10F-0 1 
0*381 30F-01 
0.48922E-01 
v. 


SLADE SURFACE 4 


THETA 

0* 73898E-01 
0.67066E-01 
0* 516 L6E-01 
0* 23609E -01 


DERIVATIVE 

-1.10561 

-1*66752 

-2.31958 

-2*82949 


2N0 DERIV 
-116.003 
- 1 12*352 
-50,7926 
*35.7093 


3 


LEAD EDGE 8-H 
TRAILING EDGE F-I 


FREESTREAH 
VELDC IT Y 
161*064 
157* 135 


MAXIMUM VALUE CRITICAL 

FOR RHQ*H VELOCITY 

24 t. 239 310.645 

241.239 310*645 


BOUNDARY A-N 
BOUNDARY &*H 


calculated PROGRAM CONSTANTS 
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PITCH HT HM1 

0.8267349E-01 0. 4 1 33675E-02 0. 1284737E-02 


HM2 

O.1353333E-02 


HM3 

0*I240588£“02 


ITMIN UMAX 

-14 25 


LAMBDA 
30* 762794 


BETA CORRECTED 
TO BOUNDARY 
48-0000 
-47,0000 



5 NUMBER OF INTERIOR ME$H POINTS = 1053 


SURFACE 


BOUNDARY values 



BV 

0 , 

1,00000 

- 0*35999 

0*64001 


blade data at intersections of vertical mesh lines with blades 


M 

0 

0*12847E“Q2 

0,25695F-02 

0*1B542E-02 

D*513S9E-02 


BLADE SURFACE 1 
TV OTOW 


0 

0*533l4E-02 
O*92485£-02 
0.127T2E-01 
0, 15945E-0 1 


0* IGOOOE 11 

3*24264 

2*B8173 

2*&Q45B 

2*31654 


BLADE SURFACE 2 
TV DTDMV 


0,B26?3I-Q1 

0»B0810E-01 

G.82671E-01 

0*B4500E“01 

G*S6313E-0| 


-D*IOOOOE 11 
1*43838 
1*42812 
1*41752 
1*40599 


STREAM sheet coordinates AND thickness table 



I* 

M 

R 

SAL 


B 


OB /DM 


1 

-Q*1L563E-Q1 

0*32335 

-0 

0. 

1 00 DOE 

-01 

*0 

24 

2 

-0* 10278E-01 

0*32385 

“0 

0* 

IOOOOE 

-01 

-0 

3 

-0*89932E-02 

0,32385 

”0 

0- 

looooe 

-0! 

-0 

4 

-0* 77084E-02 

0. 32385 

“0 

0* 

1 DODGE 

"01 

-0 


5 

-0«64237E-02 

0,37185 


0* 

L GOODE 

-01 

-0 


IM 

IV ARRAY 



I T v ARRAY 






BLADE 









SURFACE 1 

2 

3 

4 






NO* 






1 

1 


0 

19 

0000 

0000 



2 

21 


0 

19 

0000 

OQUO 



3 

41 


n 

19 

oooo 

oooo 



4 

61 


a 

19 

0000 

0000 



5 

81 


0 

19 

0000 

0000 



CO 

to 


TABLE M, - Continued. SAMPLE OUTPUT 


f H C0OR0INATFS Of INTF8SECTI0NS Of HORIZONTAL Mf$H LINES to [ FH GLADE 


MH ARRAY - 8LADE SURFACE L 


«H 

RMH 

BEH 

fl£T AH 

otowh 

0 

0.3238 

O.100OE-D1 

90,000 

0.1000E 11 

0*92251-03 

0.3238 

0 , 100OE-0T 

87.526 

3.3728 

0.2233E-02 

0,3238 

o.ioooe-ot 

43.797 

2. 9608 

O, 3 7 L3E-02 

0. 3238 

0* 1O0OE-0L 

40.472 

2.6347 

O-5394F-02 

0.3238 

o.ioooe-ol 

36*494 

2.2844 


THETA COORDINATES DF HOAIZJ'JTAL MESH LINES 

IT THETA 

-LA *0.5787 IE-01 

-13 -0.53738E-O1 

-12 -D.49604E-QL 

-il -Q,45470t-Ql 

-tO -0* A 13376-01 


IT 

ip 

FPL 

I P2 

IP i 

1P4 

A [ 1 } 

At 21 

A t 31 

A ( 4 1 

K 

H « 

i 

I T I 

* 0 








0 

i 

20 

2 

0 

2 L 

0, 

0, 

0. 

1.00000 

0.05329 

1 

2 

1 

3 

l 

22 

0, 

0. 

O, 

1,00000 

0*05329 

z 


2 

4 

2 

23 

0, 

o. 

0. 

1.00000 

0.05329 

3 

4 

3 

5 

3 

24 

0. 

0. 

0, 

1,00000 

0,05329 

4 

5 

4 

6 

4 

25 

0, 

0. 

0. 

1,00000 

0.05329 

5 

6 

5 

7 

5 

26 

0, 

0, 

0. 

1.00000 

0.05329 

6 

7 

6 

B 

6 

27 

0* 

0. 

0. 

1,00000 

0-05329 

7 

8 

7 

9 

7 

28 

0, 

0, 

0. 

1.00000 

0.05329 

a 

1 

a 

10 

8 

29 

0. 

0. 

0. 

1.00000 

0.05329 

9 

10 

9 

11 

9 

30 

0. 

0. 

0. 

u ooooo 

0.05329 

10 

11 

10 

12 

10 

31 

0. 

0. 

0, 

1.00000 

0,05329 

1 1 

12 

11 

13 

t 1 

32 

0, 

0, 

0. 
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VELOC LTV 

ANGLE ( DEG ) 




147*34 

43*69 

150.46 

40.46 

160*10 
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* 


BLADE 

SURFACE l 


* 


blade surface z 



it 

VELDCJ TY 

angle r deg i 

SURF. LENGTH 

W/WCR 

* 

VELOC LTY 

angleioegj 

SURF. LENGTH 

W/WCR 

0 

* 

0 

90-00 

0 

0 

* 

0 

-90*00 

0 

0 

0. 1235b- 02 

* 

297.81 

46.40 

0.2L52E-Q2 

0.9500 

* 

62.663 

24,90 

0*l417E-02 

0*2017 

0.2569E-02 

+ 

255.00 

43.02 

0* 3958E-02 

0*8234 


05*707 

24*02 

O.2033E-O2 

0-2759 

0 * 3854E-02 

* 

248.70 

40.15 

0.5676E-02 

0.0006 


96-491 

24*66 

0.424BE-02 

0,3106 

0.5139E-D2 

* 

245* 70 

3 7- H 

0.7321E-02 

0.7909 

* 

100*27 

24*40 

0.5660E-02 

0*3220 


1 M- 


SURFACE VELOCITIES BASED 


ON tangential components 




M 

0 

0.9225E-03 

D.2231E-0? 

Q.3713E-D2 


&LADE 
VELOCl TY 
261.23 
3 10 * 6 A 
265*50 
251.53 


surface 1 

ANGLE | DEG) 
90.00 
47,53 
43. BO 
40,47 


W/tfCR 
0.0411 
1 * 0000 
0*9549 
0.9097 



r 


BLADE SURFACE VELOCITIES 


151 


0* 


0,0100 


0.0200 


0-0300 


0-04QQ 


0.0500 


0* 

so* 

100* 

150- 200* 

250- 

300* 

350* 

400- 

450- 

500. 


1 X 

l 

l 1 

I 

+ * 

1 

1 

l 

1 


i 

l 

1 l 

1 * 

l 

L 

1 

4 

t 


l 

X l 

1 1 

1 + 

L 

1 

1 

1 

l 


1 

XI 

1 l 

+ 

l 

l 

\ 

l 

4 


1 

X 

1 1 

* 1 

1 

i 

1 

L 

1 


1 

X 

1 l 

+ 1 

1 

L 

1 

1 

1 


L 

1 

1 1 

*1 

1 

L 

1 

1 

1 


1 

IX 

1 1 

+ 1 

1 

I 

1 

1 

t 


l 

X 

1 1 

+ 1 

L 

i 

1 

l 

I 


1 

l 

1 1 

l 

l 

1 

l 

l 

1 


l 

X 

1 1 

+ 1 

1 

1 

1 

l 

1 


I 

XL 

1 1 

+ 1 

1 

l 

L 

1 

1 


1 

XL 

i 1 + 

1 

1 

1 

L 

I 

l 


X 

XI 

l 1 +. 

1 

1 

1 

L 

1 

1 


l 

l 

1 1 

t 

1 

1 

1 

1 

1 


1 

XI 

1 1 + 

L 

l 

L 

1 

1 

4 


1 

XI 

1 + L 

i 

l 

1 

£ 

1 

l 


I 

XI 

1 +1 

1 

1 

1 

1 

1 

I 


1 

X 

1 + 1 

l 

1 

£ 

1 

i 

1 


L 

IX 

l *1 

i 

l 

4 

L 

4 

1 


1 

1 X 

1 + l 

i 

l 

l 

i 

4 

1 

{ 1 

1 

1 X 

1 * 1 

i 

1 

L 

1 

1 

4 


1 

1 

l 1 

L 


l 

1 

L 

! 


1 ) 

1 X 

1 + 1 

1 

1 1 

1 

£ 

L 

1 


1 

J l X 

1 +1 

1 

1 * 

X 

1 

1 

1 

1 0 

1 

\ 1 

1 l 

1 t 

1 

1 

1 

1 

1 


1 

1 

l l 

1 

1 

1 

I 

1 

1 


1 

11 

l 1 

% 

l 

1 

1 

1 

1 


i 

J 

I 1 

1= 1 

1 

1 

1 

1 

1 


t 

1 1 

1 1 

l L 

1 

1 

1 

1 

1 


l 

L 

1 l 


1 

1 

1 

i 

1 


1 

1 1 

1 1=1 

1 

i 

£ 

1 

L 

l 


l 

1 l 

1 t 

1 

1 

I 

1 

1 

1 


1 

l t 

1 %*l 

1 

1 

1 

1 

1 

t 


1 

l ) 

X * 1 

1 

1 

1 

1 

i 

1 


1 

i r 


1 

1 

£ 

1 

1 

1 


I 

i i 

] t ( 

1 

1 

1 

1 

i 

1 


I 

i i 

1 =1 1 

1 

L 

1 

1 

i 

1 


1 

i } 

1 = 1 

L 

i 

1 

£ 

i 

1 


l 

i i 

l 1 

1 

i 

1 

L 

t 

l 


1 

i 

i 1 1 

1 

I 

1 

1 

i 

l 


I 

i 

) l L 

1 

1 

1 

1 

i 

1 


L 

L 

) 1 I 

1 

i 

1 

1 

i 

1 


l 

l 

1M i= 1 

1 

i 

1 

1 

i 

1 


I 

l 

1 l 

4 

i 

1 

L 

i 

1 

0 - 

SO* 

100. 

ISO. 200* 

Z 50* 

300. 

350 - 

400* 

450- 

500. 


VELOClfV(H) VS- MERIDIONAL STREAMLINE DISTANCE^ DOWN THE PAGE 


OO 


16 II HE * 2,9211 MJN* 
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ERROR CONDITIONS 


The error conditions are as follows: 

(1) SPLINT USED FOR EXTRAPOLATION 

EXTRAPOLATED VALUE = X.XXX 

SPLINT is normally used for interpolation, but may be used for extrapolation in some 
cases. When this occurs, the above message is printed, as well as the input and output 
of SPLINT. Calculations proceed normally after this printout. 

(2) BLCD CALL NO. XX 

M -COORDINATE IS NOT WITHIN BLADE 

This message is printed by subroutine BLCD if the M -coordinate given this subroutine as 
input is not within the bounds of the blade surface for which BLCD is called. The value 
of m and the blade- surface number are also printed when this happens. This condition 
may be caused by an error in the integer input items for the program. 

The location of the error in the main program is given by means of BLCD CALL NO. 
XX, which corresponds to locations noted by comment cards at each MHORIZ, ROOT, 
and BLCD call in the program. 

(3) ROOT CALL NO. XX 

ROOT HAS FAILED TO CONVERGE IN 1000 ITERATIONS 
This message is printed by subroutine ROOT if a root cannot be located. The input to 
ROOT is also printed. The user should thoroughly check the input to the main program. 

The location of the error in the main program is given by means of ROOT CALL NO. 
XX, which corresponds to locations noted by comment cards at each MHORIZ and ROOT 
call in the program. 

(4) DENSTY CALL NO. XX 

NER(l) = XX 

RHO*W IS X-XXXX TIMES THE MAXIMUM VALUE FOR RHO*W 
This message is printed if the value of pW at some mesh point is so large that there is 
no solution for the values of p and W. This indicates a locally supersonic condition, 
which can be eliminated by decreasing WTFL in the input. 

If RHO#W is too large, TANDEM still attempts to calculate a solution. This often 
permits an approximate solution to be obtained which is valid at all the subsonic points in 
the region. In other cases, the value of W is reduced at some of the points in question 
during later iterations, resulting in a valid final solution for these points. The program 
counts the number of times supersonic flow has been located at any point during a given 
run (NER(l)). When NER(l) = 50, the program is stopped. 

The location of the error in the main program is given by means of DENSTY CALL 
NO- XX, which corresponds to locations noted by comment cards at each DENSTY call 
in the program. 
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(5) MM, NBBI, NRSP, OR SOME SPLNO IS TOO LARGE 

If this is printed, reduce the appropriate inputs to their allotted maximum values. 

(6) WTFL IS TOO LARGE AT BLADE LEADING EDGE 

This is printed if WTFL is greater than the choking mass flow for the boundary BM. If 
this message is printed, WTFL is cut in half by the program and calculations proceed as 
usual for one outer iteration. 

(7) ONE OF THE MH ARRAYS IS TOO LARGE 

This is printed if there are more than 100 intersections of horizontal mesh lines with any 
blade surface, hi this case NBBI should be reduced, 

(8) THE NUMBER OF INTERIOR MESH POINTS EXCEEDS 2000 

This is printed if there are more than the allowable number of finite -difference grid 
points. Either MM or NBBI must be reduced. 

(9) SEARCH CANNOT FIND M IN THE MH ARRAY 

If this is printed, the value of m and the blade -surface number are also printed. The 
user should thoroughly check the input to the main program. 


PROGRAM PROCEDURE 


The program is segmented into seven main parts, the subroutines INPUT, P RE CAL, 
CGEF, SOR, SLAX, TANG, and VELOCY called by the main program TANDEM* In ad- 
dition, there are several other subroutines. All the subroutines and their relation are 
shown in figure 15* All information which must be transmitted between the seven main 
sub routines is placed in COMMON* 

Most of the subroutines in TANDEM use the same set of variables* These variables 
are all defined in the section Main Dictionary (p. 50)* All subroutines using these vari- 
ables are described prior to the main dictionary* The remaining subroutines are de- 
scribed after the main dictionary, and variables are defined with each subroutine. 

The program can handle as many as 2000 mesh points on the IBM 2-7094-7044 direct- 
coupled system with a 32 768 -word core* For 2000 mesh points to be handled an over- 
lay arrangement is used, as shown in figure 16* All subroutines not shown are in the 
main link* The total program storage requirement is 74513^ of which 46770^ is in 
COMMON blocks which are stored in the main link* The system storage requirement for 
our computer is 2764 ^ and unused storage is 300^. If there is a storage problem on 


( 8 )' 


the user's computer, the maximum number of mesh points should be reduced. The fol- 
lowing program changes are required to change the maximum number of mesh points: 

(1) Change the dimension of A, U, K, and RHO in the COMMON/AUKRHO/statement. 
This statement occurs in most subroutines. 

(2) In subroutine INPUT, change the number of values of U, K, and RHO to be ini- 
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tialized (the bound on the DO loop near statement 60). 

(3) In subroutine PRECAL, change statement 340 and format statement 1150 to re- 
flect the maximum allowable number of mesh points. Statement 340 will cause the pro- 
gram to stop if there are too many mesh points. 

(4) Change the dimensions of W, RWM, and BETA in SLAX, SLAV, TANG, VELOCY, 
and VEL. 

{5) If the number of mesh points is reduced to below 1600, the equivalence statements 
in SLAX, SLAV, TANG, VELOCY, and VEL must be changed. 

The first segment of the program is INPUT. This subroutine reads all input data 
cards, calculates constants, and initializes arrays. The next subroutine is PRECAL, 
which calculates all quantities which remain constant for a single problem. INPUT and 
PRECAL are each called once for a given problem. The remaining subroutines are 
called once for each outer iteration. The subroutine COEF calculates the entries of the 
matrix A and the vector k of equation (A7). These coefficients must be recalculated 
for each outer iteration. On the first outer iteration subroutine SOR estimates an opti- 
mum overrelaxation parameter 12 on the first call if it is not given as input. The same 
value of £1 is used for each outer iteration. SOR then finds the linear solution to equa- 
tion (A7) with fixed coefficients by successive over re taxation. Then subroutine SLAX 

calculates the streamline locations and pW„ and plots the streamline locations if de- 
ni 

sired. Subroutine TANG calculates pW Q and then pW and (S throughout the region. 
Finally, the subroutine VELOCY calculates the density p and velocity W throughout 
the region and on the blade surfaces and plots the surface velocities. 


Conventions Used in Program 

For convenience, a number of conventions are used in naming variables and assigning 
subscripts. First, several pairs of variables are spelled the same except for one letter, 
which is U in one case and L in the other. The U signifies an upper surface BC, DF, EF, 
or JK, and L signifies a lower surface ML or JL Another practice is to use the letters 
I and O in a similar manner, where I refers to the inlet or region ABMN, and O refers 
to the outlet or region FGHI. Similarly, the letter T refers to 0, and M refers to m. 
Finally, V is used to refer to vertical mesh lines, and H refers to horizontal mesh lines. 
For example, DTDMH is an array of the values of de/dm at the intersections of hori- 
zontal mesh lines with the blade. 

The variable IP is used to number all the mesh points. It starts with IP = 1 at A 
and is incremented up the vertical mesh lines one by one to the right, ending with 
IP = NIP at the last mesh point near H. The mesh spacing in the m-direction is labeled 
HM1, HM2, or HM3, and the spacing in the 8 -direction is HT. The subscript IM de- 
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notes the number of a vertical line, from IM = 1 at AN to IM = MM at GH. IT denotes 
the number of a horizontal mesh line. IT is zero along AB, increases to ITMAX at the 
highest mesh line in the region and decreases to ITMIN for the lowest mesh line in the 
region. 


Labeled COMMON Blocks 

For convenience, most variables which are used in more than one subroutine are 
placed in labeled COMMON blocks. A brief description of each labeled block is given. 
The same variable names are used in different subroutines for every variable in a 
COMMON block. The only exception is when EQUIVALENCE is used for variables in 
/AUKRHO/. The labeled COMMON blocks are as follows: 

/ENP/ is used for input variables, with the exception of those in /GEOMEN/. 

/GEOMEN/ is used for certain geometry input variables which are needed only in 
£LCD. 

/CALCON/ is used for calculated constants which are initially calculated in INPUT 
or PRECAL and do not change after this. 

/AUKRHO/ is used for the arrays A, U, K, and RHO {see section Main Dictionary) 
or the variables which are made equivalent to some of these. 

/BLCDCM/ is used for internal variables for BLCD. /BLCDCM/ is needed only to 
save certain values when overlay is used. 

/HRBAAK/ is used for variables calculated by HRB to be used in AAK. 

/RHOS/ is used to store values of p on blade surfaces. 

/SLA/ is used for streamline 6 -coordinates. 

/BOX/ is used for internal variables for the spline subroutines in order to reduce 
storage requirements. 


Subroutine INPUT 

Reading and printing of input . - The program first reads all input cards for a partic- 
ular problem. A description of the input required is given in the section Instructions for 
Preparing Input. All the input data are printed as the first output. 

Calculation of constants and initialization . - After all input is read, many of the sim- 
pler constants used throughout the program are calculated. Finally, all density arrays 
are initialized to p| n (RHOEP). 
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Subroutine P REGAL 


Calculation of constants. - The prerotation X and the average relative velocity W le 
at the blade leading edge are calculated first by an iterative process (eqs. (B7) to (B9)). 
During this calculation the input weight flow (WTFL) is checked to see if it is larger than 
the upstream choked flow value. IE so, it is cut in half and the computation of X and 

* s re P ea t^d. Maximum values of the mass flow parameter pW (eqs. (BIO) to (JB 12)) 
and the critical velocity W cr are then calculated at the leading and trailing edges of the 
blade row. The flow angles and /3 Qut are also computed at the upstream and down- 
stream boundaries (appendix B) from the values, /3 lg and 0 te , given at the leading and 
trailing edges. 

Calculation of vert ical m esh line arrays . - BLCD is called for the four blade sur- 
faces obtaining 8 -coordinates (TV) and slopes (DTDMV) where the vertical grid lines 
meet the blades. By using the TV array, the integer arrays (ITV and IV) are calculated. 
Finally, by using DTDMV, the blade-surface angles (BETAV) are calculated. 

Calculation of horizontal mesh line arrays. - MHORIZ is called once for each of the 
four blade surfaces to obtain the m -coordinates (MH) and slopes (DTDMH) where horizon- 
tal grid lines meet the blade surfaces. Then by using cubic spline interpolation (SPLINT) 
and the MH array, RMH and BEH are calculated. Finally, the blade-surface angles 
BETAH are calculated by using DTDMH. 


Subroutine COEF 

Subroutine COEF controls the calculation of the finite -difference coefficients of u 
in equations (A2) to (A6) (elements of the matrix A in eq. (A7)). At the same time, it 
computes the constants of the finite -difference equations (components of k in eq. (A7)). 

Calculating coefficients and constants thro ughout region . - COEF progresses from 
left to right through the blades. COEFP and COEFBB are called along each vertical 
mesh line for the calculation of the coefficients and constants. COEFP is called in the 
periodic regions upstream and downstream of the blades, and between front and rear 
blades for the nonoverlapping case. COEFBB is called in the regions between upper and 
lower blade surfaces. 

Corrections to co efficients and constants . - At certain points in the solution region, 
corrections must be made to the coefficients and constants calculated by COEFP and 
COEFBB. This is done at the end of COEF if points B, J, C, E, or F (see fig. 4) on the 
blade surfaces coincide with mesh points in the solution region. Corrections are also 
made along line KL and the line to the right of CD. The periodic boundary condition equa- 
tions are applied herein (see eqs. (A5) and (A6) and the explanation following them). 
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Subroutine COEFBB 


Subroutine COEFBB computes the coefficients and constants k. along a vertical 
mesh line from blade to blade. It has a second entry point (COEFP) with completely sep- 
arate code, which computes coefficients and constants in periodic regions. Both COEFP 
and COEFBB proceed up a vertical mesh line, one point at a time. 

In both the periodic and the blade -to-b lade cases, HRB is called initially to compute 
the values of h, r, and b required in equation (A2). These values are then altered for 
special cases. In COEFBB they are altered along lines CD and KL, and in COEFP along 
periodic boundaries. COEFBB also calls BDRY12 and BDRY34 to obtain special values 
of h, r, and b when mesh points are within one mesh space of a blade boundary. Fi- 
nally, both COEFP and COEFBB call AAK to compute A and k from equations (A2). 


Subroutine HRB 

Subroutine HRB calculates values of h, r, and b for use in equations (A2). Each 
time it is called, it computes these values for a single mesh point. 


Subroutine AAK 

Subroutine AAK is called by COEFBB and computes the coefficients and the con- 
stants kj of equation (A2) at a single point. 


Subroutine BDRY12 

Subroutine BDRY12 is called by COEFBB. It alters the values of h and r calcula- 
ted by HRB for point 1 or 2 (see fig. 17) if either of these points lies on a blade surface. 
It also defines the constants KAK and KA used to alter A and k in COEFBB or COEF. 


Subroutine BDRY34 

Subroutine BDRY34 is called by COEFBB. It alters the values of h, r, and b cal- 
culated by HRB for point 3 or 4 (see fig. 17) if either of these points lies on a blade sur- 
face. It also defines the constants KAK and KA used to alter A and k in COEFBB or 
COEF. 
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Subroutine SOR 


This subroutine solves the finite-difference matrix equation (A7) by the method of 
successive over relaxation (ref. 10). The same section of code is used both for calcula- 
ting the optimum overrelaxation factor Q and to solve equation (A7). If a value of ORF 
greater than 1 and less than 2 is given as input, it is used for the overrelaxation factor. 
Otherwise a value is estimated by the program. 

In equation (A8), the subscript i denotes the index of an unknown mesh point. In the 
program, i is replaced by IP. The subscript j in equation (A8) denotes the index of 
neighboring unknown mesh points. For each i, there are only four values of j for which 
a^. is nonzero, which are the negative values of the coefficients A(IP,1), A(IP, 2), 

A (IP, 3), and A (IP, 4). The value of j is determined by the index of the proper neighbor- 
ing point. These indexes are named IP!, IP2, IP3, and IP4. These indexes are defined 
so that has the coefficient A (IP, 1); the other indexes are defined similarly. 

Estimation of optimum overrelaxation factor . - H ORF = 0 as input, the optimum 
value for the over re taxation factor SI is estimated on the first outer iteration by using 
equations (B3) and (Bl) of reference 11. Equation (A8) is used to calculate u m+ ^ from 
u m for equation (B3) of reference 11, with 0=1, and k = 0. Equation (A8) becomes 


u m+1 

l 


i-1 


-Evr-E 


J=1 


j=i+ 1 


a. u m 
i] J 


( 6 ) 


To start, u? = 1 for all i. The maximum value of the ratio u™ + ^ / u™ is calculated for 
a given m and is given the name LMAX. After convergence, the optimum value of the 
overrelaxation factor O can be calculated by O = 2/(1 + -^1 - LMAx). This procedure 
is explained in appendix B of reference 11. 

Solution of matrix equation by subrouti ne SOR . - With a value of O either as input 
or estimated by the program, equation (A8) can be used iteratively to calculate a se- 
quence {u m j that will converge to a solution to equation (A7). During each iteration, 
the maximum change of the stream function is calculated. When this maximum change is 
reduced below 10 , the iteration is stopped, and the current estimate of the stream 

function is accepted as the solution. 


Subroutine SLAX 

Subroutine SLAX, by calling subroutines SLAVP and SLAVBB (entry points of SLAV), 
computes the meridional mass flow component pW m at all points on vertical mesh lines. 
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Subroutine SLAX also calculates and plots streamline locations. 

Calculating pW m throughout region. - Subroutine SLAX progresses from left to 
right"through the blades, calling SLA VP and SLAVBB along each vertical mesh line. 

SLA VP is called in the periodic regions upstream and downstream of the blade and be- 
tween blades for the nonoverlapping case. SLAVBB is called in the blade -to -blade 
regions. 

Plotting streamlines . - When subroutine SLAX reaches the right end of the region, 
all information is available from SLA VP and SLAVBB for the streamline plot. The plot- 
ting printout is done by PLOTMY, which, with the necessary further subroutines PISTUG 
and KHAR, is described completely in reference 12. 


Subroutine SLAV 

Subroutine SLAV has two entry points, SLA VP ana SLAVBB. SLAVP is called in 
’periodic regions, SLAVBB from blade to blade. Both entry points make use of a common 
section of code at the end of SLAV. 

Calc ul ation of 5u/9<3 a nd st ream line locations. - SLAVP and SLAVBB compute 
3u/36? along each vertical mesh line. The derivative du/d9 is estimated at each mesh 
point from a cubic spline curve (SPLINE) of the stream function u. 

SLAVP and SLAVBB also calculate values of 9 corresponding to given values of the 
stream function. These values are printed out and are also used for the streamline plot. 
The stream function is a one-to-one function of distance in the 6 -direction along most 
vertical mesh lines. Therefore, cubic spline interpolation (SPLINT) can be used to ob- 
tain 9 as a function of u. 

Calculation of pW m - - SLAVP and SLAVBB use the derivatives 3u/S0 to calculate 
P W m at each mesh point. The equation 3u/8(? = brpW m /w (eq. (3)) is used. Values of 
pW m are stored in RWM for interior mesh points, and in WMB where the blade surfaces 
are intersected by vertical mesh lines. 

Calculation of mass flow parameter pW on blade surfaces . - Where each vertical 
mesh line meets a blade surface, pW is calculated from PW m by the equation 



Subroutine TANG 

Subroutine TANG calculates the tangential mass flow component pW g at all points 
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on horizontal mesh lines. This process is complicated by the fact that the horizontal 
mesh lines are shifted in crossing the boundary KL. 

Location of points on horizontal mesh lines . - Subroutine TANG begins at the bottom 
line of the region and proceeds upward to the top of the region, moving from left to right 
along each horizontal mesh line. On a given mesh line, the first point in the region is 
located by comparing IT for that mesh line with ITV for each of the four blade surfaces of 
successive points along the line. After an initial point in the region is located, TANG 
moves to the right along the line until it encounters the downstream boundary GH or a 
blade surface. Once again, TANG locates a blade boundary by comparing IT with ITV of 
the blade surfaces. ROOT is called to calculate mesh spacing at the end points when they 
are located on one of the blades. TANG stores the meridional distance and stream - 
function value of each point located along a line into the arrays SPM and USP. 

Calculation of pWg. - When a horizontal mesh line exits from the solution region, 
subroutine SPLINE is called with SPM and USP to calculate 9u/Sm at each point along 
the line. The product pW^ is then calculated from 9u/9m by using 9u/9m = -bpWg/w 
(eq. (2)). 

Calculation of mass flow parameter pW_and flow angles at interior points . - At each 
interior point, pW is calculated by 

pw ^/(pw m ) 2 + (pW fl ) 2 
and the angle 0 is calculated by 


tan /3 


PW, 


pw 


m 


These values are stored in W and BETA for all interior points. 

Calculation of mass flow par ameter pW on blade surfaces . - Where each horizontal 
mesh line meets a blade surface, pW is calculated from pW Q by the equation 
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Subroutine SEARCH 

Subroutine SEARCH is used by TANG in the calculation of the mass flow parameter 
pW on blade surfaces. The distance (DIST) corresponds to some element in the MH ar- 
ray for a particular surface. SEARCH locates that element and returns its subscript to 
TANG. TANG then uses a corresponding element in the BEH array in calculating pW. 


Subroutine VELOCY 

Subroutine VELOCY calculates densities p and velocities W from the mass flow 
parameter pW at all points throughout the solution region and on the blade surfaces. It 
also plots the surface velocities. 

Solving for densities an d velocities throu ghout region . - VELOCY progresses from 
left to right through the blades, calling VELP and VELBB for each vertical mesh line. 
VELP is called in the periodic regions, and VELBB is called from blade to blade. When 
the right boundary of the solution region is reached, VELSUR is called once to compute 
the blade-surface velocities. 

Plotting of velocitie s. - After VELOCY calls VELSUR, all information is available 
for the plot of surface velocities. The velocities are plotted by using different symbols 
for front and rear blades, upper and lower surfaces, and velocities based on both merid- 
ional and tangential components. Velocities based on meridional components are plotted 
if j (3 J ^ 60°, and velocities based on tangential components are plotted if J {3 J s 30°. 
Plotting is done by PLOTMY, which is described in reference 12. 


Subroutine VEL 

Subroutine VEL has three independent entry points, VELP, VELBB, and VELSUR. 
VELP and VELBB compute velocities in the periodic and blade -to -blade regions, and 
VELSUR computes velocities on the blade surfaces. None of these entry points share 
common code in VEL. 

The maximum relative change in density p along a blade surface is calculated in 
VELBB and VELSUR and is called RELER. If RELER is less than 0. 001, the outer iter- 
ation is considered to be converged, and the calculations are stopped on the following 
iteration. 

Calculation of p and W. - Both VELP and VELBB proceed from left to right through 
a region, and upward at each vertical mesh line from boundary to boundary. VELSUR 
proceeds along the four blade surfaces one at a time. VELP, VELBB, and VELSUR cal- 
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culate density and velocity from pW by calling the DENSTY subroutine at each mesh 
point and boundary point. Along the blade surfaces, VELBB and VELSUR also calculate 
the ratio W/W cr . 

Printing of velocities . - VELP and VELBB print interior velocities and flow angles 
as they are calculated. Surface velocities, blade- surface angles, arc lengths, and ratios 
of velocity to critical velocity are printed at the end of VELSUR. 


Subroutine BLCD 

Subroutine BLCD calculates the 6 -coordinate and de/dm of a blade surface for any 
given value of m. There are four entry points to BLCD corresponding to the four blade 
surfaces. 

The first time that BLCD is called for a particular blade surface, the coordinates of 
the first and last spline points are calculated. These points are tangent to the leading - 
and trailing-edge radii, respectively. The parameters defining the spline curve are also 
calculated at this time. 

Each blade surface is defined by the leading- and trailing-edge radii and by a cubic 
spline curve, which is a piecewise cubic polynomial. The procedure is to scan the spline 
points to determine which interval the m -coordinate is in and then to calculate the 0- 
coordinate and derivative. 

The arguments for the entry points of BLCD are defined so as to be called by ROOT 
to determine the m -coordinate of an intersection of a horizontal mesh line with the blade. 
Most of the information needed by BLCD is in labeled COMMON blocks. These variables 
are found in the main dictionary. 

The input argument is 

M meridional streamline coordinate, m 

The output arguments are as follow: 

THETA 9 -coordinate of blade surface at m 
DTDM d9/dm of blade surface at m 

INF used when d(?/dm is infinite; INF is normally 0, but set equal to 1 if d0/dm 

is infinite 
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Function IPF 


A mesh point in the solution region can be numbered in one of two ways. The first is 
by coordinates of mesh line intersection, IM and IT. IM is the number of the vertical 
mesh line, beginning with 1 at the upstream boundary AN. IT is the number of the hori- 
zontal mesh line, beginning with 0 at the leading edge of the upper surface of the front 
blade. The second numbering system is by point count, using EP. IP increases up each 
succeeding vertical mesh line from left to right through the solution region. IPF returns 
the value of IP corresponding to given coordinates, IM and IT. 


Main Dictionary 


The Main Dictionary applies to ail the previously discussed subroutines. 


i 

A 

A12, A34 

AA 

AAA 

AANDK 

AATEMP 

ADD 


ADDL 

ANS 

AR 

AZ 

B 

B12, B34 
BB 


array of coefficients of u (i. e. , elements of a- of matrix A in 
eq- <A7)) 

a-jg, a^ in eq. (A2) 

temporary variable in PRECAL and BLCD 
array used for temporary storage 
see Input 

temporary location for AANDK in SOR 

logical variable in TANG, indicating need to add 1 to stream func- 
tion at a mesh point prior to spline fit of stream function along 
a horizontal mesh line 

logical variable in TANG, indicating entrance into region where 
ADD applies 

result of calls on ROOT in TANG and DENSTY in VEL 
see Input 
aQ in eq. (A2) 

array containing stream -channel thickness b at the four points 
adjacent to a point for which AAK is called 

b 12’ b 34 111 eq ' < A2 > 

temporary variable in PRECAL and BLCD 
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BE 

BEH 

BESP 

BETA 

BETAH (BETAV) 

BETAI (BETAO) 
BETI (BETO) 

BLDAT 

STAIN 

BTAOUT 

BV 

BZ 

CDMBIT (CDMBOT) 
CHANGE 

CHORD 

CMM 

CP 

CPTIP 

DBDM 

DELTV 

DIST 


array of values of b 

array of values of b 
blade surfaces 


at vertical mesh lines 

where horizontal mesh lines meet the four 


see Input 

array of values of j3 at interior mesh points 

array of values of 0 where horizontal (vertical) mesh lines meet 
the four blade surfaces 


see Input 

array of angles at tangent points of leading- (trailing-) edge radii 
with the four blade surfaces (see input BETH, 2, 3, 4 and 
BETOl, 2, 3, 4) 

see Input 

free -stream angle at upstream boundary AN based upon /3j g , 

calculated by eq. (B14) 

free -stream angle 0 ^ at downstream boundary GH based upon 

calculated by eq. (B14) 

array of stream -function boundary values on the four blade sur- 
faces 


stream -channel thic kn ess at a point for which AAK is called 

temporary grid locations along meridional axis in INPUT 

change in value of stream function at a particular point during an 
iteration of SOR 

array containing the meridional chord distances of each of the 
four blade surfaces (see input CHORDF and CHORDR) 

temporary variable in BLCD 

specific heat at constant pressure, c p 

^ c p^in 

array of slopes at vertical mesh lines of spline curve for stream - 
channel thickness 


increment in Q -coordinate in VEL 

meridional distance in SEARCH from a blade leading edge to 
where a horizontal mesh line meets a blade surface 
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I 


DMLR 


DTDM 

DTDMH (DTDMV) 

DTLR 

DUDM 

DUDT 

EM 

EMK, EMKM1 
ERROR 

ERSOR 
EX PON 
FIRST 
GAM 
H 

HM1 

HM2 

HM3 

HT 

I 

11,12 


tolerance for mesh points near a boundary in m- direction (If a 
mesh point is closer than DMLR to a boundary, the point is 
considered to be on the boundary. ) 

de/dm along a blade surface in BLCD 

array of d0/dm where horizontal (vertical) mesh lines meet the 
four blade surfaces 

tolerance in 9 -direction (see DMLR) 

array of derivatives of stream function du/dm along horizontal 
mesh lines in meridional direction 

array of derivatives of stream function du/d£? along vertical mesh 
lines in 6 -direction 

array of second derivatives of spline curves for each blade sur- 
face, calculated by SPLN22 in BLCD 

temporary variables for EM in BLCD 

maximum absolute value of change in u at any point for an over- 
relaxation (SOR) iteration 

see Input 

i/(y - l) 

initial value of some index 

see Input 

array containing mesh spacing h between the point for which 
AAK is called and the four points adjacent to it 

mesh spacing in m -direction from upstream boundary through 
front blade 

mesh spacing in m-direction for overlapping portion of front and 
rear blades, or between blades for the nonoverlapping case 

mesh spacing in m-direction through rear blade to downstream 
boundary 

mesh spacing in 8 -direction from blade to blade 

temporary integer variable in PRECAL, SLAX, SLAV, and 
SEARCH 

temporary integers in SLAV 
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IEND 


IH 


IHS 

IM 

IM1 (EMT) 

IM2 

IMS 

IMSL 

IMSS 

IMTM1 

INF 

INIT 

ENTU 

INTVL 

IP 

IP1, EP2, IP3, IP 4 

IPCD, IPKL 
IPL (IPU) 

IPLM1 {IPUP1) 

IS 


integer variable set equal to I when final convergence to a solution 
is reached in the outer iterations on a given set of data 

array containing current number of intersections of horizontal 
mesh lines with each of the four blade surfaces as intersections 
are located 

integer variable in BDRY34 and TANG for counting intersections 
of horizontal mesh lines with blade surfaces 

index of mesh line in meridional direction (ra -direction) 

integer variable in TANG indicating the vertical mesh line index 
of the first (final) point in the region of a horizontal mesh line 

IM1 + 1 

array containing total number of intersections of horizontal mesh 
lines with each of the four blade surfaces 

temporary variable in PRECAL 

temporary variable in PRECAL, VELOCY, and VEL 

IMT - 1 

variable in PRECAL indicating when an infinite slope is located at 
a blade leading- or trailing- edge in a call on BLCD 

array used to indicate whether BLCD has been called previously 
on a given blade surface 

temporary integer streamline value in SLAV 

see Input 

index of iffesh point 

value of IP at the four adjacent points to the mesh point under 
consideration 

temporary IP along lines CD and KL in COEF 

value of IP where a vertical mesh line meets a lower (upper) sur- 
face or boundary 

value of IP on a vertical mesh line adjacent to a lower (upper) 
surface in VEL 

integer variable in SEARCH for indicating where a horizontal mesh 
line intersects a blade surface 
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IT 

IT3, IT4 

ITER 

ITI 

ITMAX (ITMIN) 

ITO 

ITV 


IT VI, ITV2 
IT VIM 1 

itvl (itvu) 

ITVLP1 
ITVM1 (ITVP1) 

ITVUM1 

IV 

IVMM 

J 

K 

KA 

KAK 

KK 

KKK 

L 

LAMBDA 


index of mesh line in $ -direction 

value of IT for the adjacent points (3 and 4) to mesh point under 
consideration 

outer iteration counter 

horizontal mesh line index in TANG one period below IT, IT-NBBI 

maximum (minimum) value of IT in mesh region 

value of IT at origin of coordinates at leading edge of front blade 

array of horizontal mesh line indexes (IT) corresponding to inter- 
sections of vertical mesh lines with blade surfaces (ITV{IM, 
SURF) is the IT value for the mesh point in the region on verti- 
cal mesh line EM which is closest to blade surface (SURF). ) 

temporary storage of ITV in TANG 

temporary ITV in TANG 

ITV of the lower (upper) blade surface on a given vertical mesh 
line 

ITVL + 1 

ITV of a blade surface in COEFBB for the vertical mesh line to 
left (right) of line under consideration 

ITVU - 1 

array containing value of IP at the base of each vertical mesh line 
temporary storage of IV in COEF 

temporary integer variable in INgUT, SLAX, and SLAV 
array of constants; vector k in eq. (A7) 

integer array indicating which of the four points surrounding a 
mesh point lie on a boundary 

real array giving the stream -function values of boundary points 
surrounding a mesh point 

integer counter in BLCD 

array containing information used in plotting subroutine PLOTMY 
temporary integer variable in SLAV 
X 


54 



LAST 

LER 

LMAX 

LOC 

LOWER 

M 

MB I 

MB 12 

MBI2M1 

MBJ2P1 

MBXM1 

MB IP X 

MBIT, MBOT 

MBO 

MB02 

MB02M1 

MB02P1 

MBOM1 

MBOP1 

MH 

MLE 

MM 

MMLE 

MMMI 

MMMSP 

MR 


final value of some index 

array indicating location of error messages printed toy program 

maximum value of u™ + ^ j u P for eq. (B2) of ref. 11 

integer variable in SLAV specifying which entry point (SLA VP or 
SLAVBB) was used 

integer variable representing one of the lower blade surfaces, 

2 or 4 

meridional coordinate, m 

see Input 

see Input 

MB 12 - 1 

MBI2 + 1 

MBI - 1 

MBI+ 1 

temporary grid locations along meridional axis 

see Input 

see Input 

MB02 - 1 

MB02 + 1 

MBO - 1 

MBO + 1 

array of m -coordinates of intersections of horizontal mesh lines 
with the four blade surfaces 

array of m -coordinates of leading edges of the four blade surfaces 
(see input MLER) 

see Input 

temporary meridional distance in BLCD 
MM - 1 

temporary meridional distance in BLCD 
see Input 
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MRTS 

MSL 

MSP 

MSPMM 

MV 

MVIMl 

NBBI 

NBL 

NER 

NI 

NIP 

NP1, NP2 

NRSP 

NSP 

NSPI 

NSP Ml 

OMEGA 

ORF 

ORFOPT 

ORFTEM 

P 

PITCH 

R 


integer switch in PRECAL indicating when infinite derivatives 
would be encountered in a call on MHORIZ 

temporary storage for MV array during plotting in SLAX 

array of m -coordinates of spline points for each blade surface 
measured from its leading edge (see input MSPl, 2, 3, 4) 

temporary meridional distance in BLCD 

array of m- coordinates of vertical mesh lines 

temporary value of MV in TANG 

see Input 

number of blades 

array indicating number of times certain error messages are 
printed by program 

number of streamlines blade to blade in SLAV 
number of interior mesh points 

integer counters in VELOCY indicating number of plotted blade- 
surface velocities 

see Input 

number of spline points 

array containing number of spline points on each of the four blade 
surfaces (see input SPLN01,2, 3, 4) 

NSP - 1 

see Input 

see Input 

upper bound for estimate of optimum from eqs. (Bl) and (B2) 
of ref. 11 

temporary storage for ORFOPT 

array containing information used in the plotting subroutine 
PLOTMY 

2jt/NBL, 8 -coordinate from blade to blade 

array of densities p at the four points adjacent to a point for 
which AAK is called 
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RATIO 
RE LEE 

RHO 

RHOl 

RHOB 

RHOHB 

RHOEP 

RHOMBI 

RHOMB 2 

RHOMM 

RHOT 

RHOVB 

RHOVI 

RHOVO 

RHOWMI 
RHOWMO 
RI (RO) 

RM 

RMDTL2 (RMDTU2) 
RMH 

RMI (RMO) 


value of u™ + for use in eqs. (B2) and (B3) of ref. 11 

maximum relative change in density at surface mesh points, be- 
tween two outer iterations 

array of densities p at interior mesh points 

average density p at upstream boundary AN 

temporary storage in VEL for a value of p on a blade surface 

array of densities p at horizontal mesh line intersections with 
the four blade surfaces 

see Input 

average density p at leading edge of front blade 
average density p at trailing edge of rear blade 
average density p at downstream boundary GH 
temporary value of density p 

array of densities p at vertical mesh line intersections with the 
four blade surfaces 

average value of pW at front-blade leading edge or upstream 
boundary AN 

average value of pW at rear-blade trailing edge or downstream 
boundary GH 

maximum value of pW at leading edge of front blade 

maximum value of pW at trailing edge of rear blade 

array of leading- (trailing-) edge radii on the four blade surfaces 
(see input RI1, 2 , 3, 4 and R01,2, 3, 4) 

array of r -coordinates of the mean stream surface radii at verti- 
cal mesh lines 

(r de/dm) at vertical mesh line intersections on lower (upper) 
blade surfaces 

array of r -coordinates of the mean stream surface radii where 
horizontal mesh lines meet the four blade surfaces 

array of r -coordinates of mean stream surface radii at the inlet 
(outlet) of the four blade surfaces 


RMM 


temporary meridional distance in BLCD 


RMSP 

see Input 

HWM 

array of pW m where vertical mesh lines intersect the four blade 
surfaces 

KWT 

array of pWg where horizontal mesh lines intersect the four 
blade surfaces 

RZ 

density Pg at point for which AAK is called 

S 

meridional distance between two adjacent blade- surface spline 
points in BLCD 

SI (ST) 

blade -surface number at beginning (end) of a horizontal mesh line 
in TANG 

SAL 

array of values of sin a = dr /dm at each vertical mesh line 

SIGN 

integer constant in BLCD 

SLCRD 

see Input 

SPLNO 

number of input spline points on a blade surface 

SPM 

array of m -coordinates along a horizontal mesh line in TANG 

SRW 

integer code variable that will cause certain subroutines to write 
out useful data for debugging: 

If SRW =13, SPLINE will write input and output data. 

If SRW = 16, SPLINT will write input and output data. 

B SRW = 18, SPLN22 will write input and output data. 

E SRW = 21, ROOT will write input and successive estimates of 

the root to which it is converging. 

STGR 

array of 8 -coordinates of center of each trailing -edge radius with 
respect to the center of its leading -edge radius (see input STGRF 
and STGRR) 

STRFN 

see Input 

SURF, SURFBV 

integer variables referring to one of the four blade surfaces 

SURFL 

array of blade -surface lengths at vertical mesh line intersections 
for each of the four blade surfaces 

SURVL 

see Input 

Tl,T2 

elapsed time in clock pulses (1/60 sec) 
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TBI 

tan %e 

TBI1, TBIT 

temporary TBI 

TBO 

tan ^te 

TBOM, TBOT 

temporary TBO 

TGROG 

2 yR/(y + 1) 

TH 

6 -coordinate from leading edge of front blade to a horizontal mesh 
line 

THETA 

9 -coordinate of a point along a blade surface in BLCD 

THK, THKM1 

temporary variables in BLCD 

THLE 

array of 6 -coordinates from origin of front blade to leading edge 
of each blade surface (see input THLER) 

THSP 

array of 9 -coordinates of spline points for each blade surface 
measured from its leading edge (see input, THSP1, 2 , 3, 4) 

TIME 

elapsed time in minutes 

TINT 

array of 9 -coordinates in SLAV where plotted streamlines cross 
vertical mesh lines 

TIP 

see Input 

TPP 

T" 

TSL 

array of 0 -coordinates of plotted streamlines 

TSP 

array of 8 -coordinates of points along a vertical mesh line in 
SLAV 

TTIP 

T/Tl 
' in 

TV 

array of 9 -coordinates where vertical mesh lines meet the four 
blade surfaces 

TWL 

2u)X 

TWLMR 

2uA - (wr)^ 

TWW 

2w/w 

U 

array of stream -function values at each mesh point, or of the 
eigenvector associated with spectral radius p(Lj), as estimated 
by the power method (ref. 11) 

UINT 

array of values of stream function for which it is desired to obtain 
interpolated values of 6 -coordinate in SLAV 
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UNEW 


UPPER, UPPRBV 
USP 
VI (VO) 
w 

WCR 

WCR1 (WCRO) 
WMB 

WTB 

WTFL 

WTFLSP 

WWCRM 

WWCRT 

XDOWN 

YACROS 


new value of stream -function estimate at a single point, calculated 
by eq. (6) 

integer variables representing one of the upper blade surfaces, 

1 or 3 

array of values of stream function along a vertical or horizontal 
mesh line, including boundary points 

average relative velocity at the leading (trailing) edge of the front 
(rear) blade 

array of relative velocities W at unknown mesh points, also used 
for storing pW 

critical velocity on a blade surface 

critical velocity at leading (trailing) edge of front (rear) blade 

array of pW m where vertical mesh lines intersect the four blade 
surfaces 

array of pW^ where horizontal mesh lines intersect the four 
blade surfaces 

see Input 

see Input 

array of ratio of blade -surface velocity (based on meridional com- 
ponent) to critical velocity 

array of ratio of blade- surface velocity (based on tangential com- 
ponent) to critical velocity 

array of m -coordinates where surface velocities are plotted 

array of surface velocities to be plotted 
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Program Listing for Subroutines Using Main Dictionary 


COMMON ,NEftUl 

COMMON /AUiUHJ/ A i 2000 ,4* , U( 2000 I , K 12 OOO ) , RHG< 23 3 3 i 
CtiMMUN /INP/ GAM ,AR , TIP ,RNUlP,ttTFL ,MTfLSP , OMEGA, ORF, BETA I f BET AQ, 
LM BI, M 80, M B 12 ,MBQ 2, rtf* , N8Bi ,NB L ,NR 5P , MR I 50 i ,RMSP153J ,8ESPf53 >, 
2BL0AT . AAXOR, ER SUR , STRF N , SLCftD , 1 H TV 1. , SURVL 
CJrtMJN /CALCUM/ M3lHl,M3IPI,M3 0Ml,rtB0Pl,MBIZrtl , MSI 2 PL , MBQ2M1, 
IMBDZPL.MMMI, HM L t HM2 iHH 3,HT ,£>TLR ,DMLR ,PI TCH ,C R, EXPQN ,TrfW , CPT 1 P f 
2TGKQG, TUlrTBG,LAMODA,TWL f I THIN*! TMA X , All P tl MS 1% } , BY * A ) f MV i 103 * t 
31 VI lul It ilVC 10 0* 4* ,TVI 100,4 J tDTDMVtlOO ,41 , BET AV <13 0 *4 J , 

4MH< 1 0 0 ■ 4 J * O T DM Hi t 100 , 4 J ,Bt TAH (1 00 *41 ,ftMHU30,4l ,BEH<103,4|, 

SRHI lUQI,3EU00) tDBDMUOOI ,$ALUO0» .AAAC100) 

COMMON /GEUMIN/ CHJRJ44I ,S?bRI4l ,HLE 141 , T H LE < 4 I ,RNI (4* ,RMQI4J, 
1KU 4j f RJt 41, BETI I4J ,5£TUI41 ,NSPI 14) , MSP (58 ,4) ,THSP<S3 ,4) 

COMMON /RHOS/ RHUHa<lOO,4*,RrtOVBll0Q,4) 

COMMON /ocCuCM/ £Mt50,4t > !NE T(4l 

INTEGER aLDAT,AANDK,ERSDR*$TRf N»$LCRO, SURVL, A ATEMP, SURF, SURFBV, 
LFIRST, UP? fcK, UPPRbV, Si, ST, SRW 

REAL R,kAR»L AMBDA,LrtAX,Mri, MLE , MR ,H5L ,HSP , MV , MV I Ml 
CALL TIUEllTlj 
10 ItHQ = -i 
ITER ^ 0 
DU 20 SU«F = 1,4 
20 1 M i T i SURE I = 0 
CALL INPUT 
CALL P REGAL 
30 CALL CUEF 
CALL 5 UR 
CALL TirtEHTZ* 

T i rt L s 4T2-in/3oiH>. 

HR IT fc f 6f LOGO ) TIME 
CALL SLAX 
CALL TANG 
CALL VLL3CV 
CALL TIMEHT2) 

T IM (T^-Tl J/36O0, 

NRITfcl 6, 1000 » TIME 
IF<NER<2}*GT*0J GO TO 10 
IP l £ END I 30,50,10 

100 0 FORMAT (BHiUME = ,F7.4,5H MI N» \ 

mo 


SUBROUTINE INPUT 

INPUT REALS AND PRINTS ALL INPUT DATA CARDS AND CALCULATES HUR 12 ONT AL 
SPACING <MV ARRAY I 

CJMMJN SRrt,l TER, I END ,LER 12 1 ,N£R<2) 

COMMON /AUKLHU/ Af2Du0,4l ,0120001 ,K<2003I ,RrtOl23331 
COMMjN /IMP/ GAM, AR f T IP , R rfUl P , H TFL , WTF LSP , OMEGA, ORF , BET A l , SET AO, 
IMS I, M BO, *B I2,rtB02,MM»N8Bl , HB L , NR $P ,MR { 50 1 ,RMSP (50) ,B£S PI S3 I, 
2SL0AT, AANDX, kR SUR, S RF N , SLCRD , I N TV L , SURVL 
COMMON /GALLON/ Mb 1 M i , MEH P 1 , MB UMl ■ MB0P1 , Mb 12 Ml , MB! 2 PI t MB02HI , 
IMROZPlfMMMl, HMl,HM2,riH3 f riT,0TlR t DMLR,PITCH,CP,EXPUNprriW,CPTlP f 
2TGRDG, T8I ,T8U, LAMBDA, TNL, I TWIN, I TMAXtNiP,! MS 14 1 , BV I 4 i , MV t 103 * , 

31 VC 1011, ITVI 100,41 , IVC 100, 4t ,OT0MVtiOO,4i , BET AV (13 0 ,4 J , 

4MH( 1 U 0 , 4 1 , 0 f DMH I 100,41 . 3 E TArt ( 100 ,4 i ,KMHUO0,4J ,8£Na33 f 4», 

5Rrt( lUOl,3fetlOO|,ua0MClOOI ,SAU100I ,AAA1L00I 
COMMON /GEQM IN/ CHORD t 4J ,STGH(41 ,MLfc 14 i ,THLE 14 J ,RMI(4l ,RM0I4|, 

1RII 4 1 , KU l 4>,BETH4i,5ETiiI4l ,N$PI (4 i , MSP (SO ,4) ,THSP{S3,4l 
CUMMUN /RHUS/ R HUNtH 100, 4 1 ,RHO VB U 00 I 

INTEGER 8L DA F, A AND*v, ER SUR , STRF N , SL€ RO , SURV L , AATEMP ,SU RP , SO RF8V , 
1F1RST, UPPER, UPPRBV,Sl , ST,SRW 
REAL K ,RAR,L AMBOA , LM A X ,MH , MLE,«R ,MSL , MSP , MV , MV I Ml 
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READ A HO PRINT ALL INPUT 0 A TA 


WRITE! 6, 1000) 

KEAU ( 5, I LOO 1 
WRITE! 6, 110Q I 
WRITE! t, lllQ * 

READ !5, 1030) GA M , A ft , T I P , Rn U ! P , W TF L # It F F LS P , OH£ G A , QKF 
WRITE! 1Q4U I GAM,AK,TIP,ftHUlP .WTFL , WTFLS P , OMEGA , QRF 
WRITE! 6, 1120 I 

READ I 5f 103UiStTA l tB£ TAO , CHORD 1 1 1 ? STbR 11 1 ,CHQROI3t ,$T GR( 3 ) * 
!HLtl3)tTHLE| 3* 

*ft | Tfc! fe, 10 AO IB ETA I, BE TAU , C HUKU t 1 J ,STGR 11 > ,CHURD!3) ,$TGR!3 ), 
1HL t ( 31 ,THLt! 31 
WRITE! 6, i 130 I 

READ I 5,1010) MBI »MBJ ,MBi2,MB02 ,MMtNBBl ,NBL,Nft$P 
WRITE! 6, I GIG 1 NBi ,MBO,MbI 2 v HB02 , MM , NBbJ « NBL , MRS P 
DO 10 J = 1 f 4 

IF I J * £Q . i 1 WftITE!6,ll4Q) 

IF ! J .EU.21 Wft I TE £ 6, 1 1 50 1 
IF !J.£U,3) Wft HE ! 6,iib0) 

IF U.EU.4I WK 1 TE I 6, 1 1 70) 


WRITE! 6,1 1 BO ) JfJtJiJfJ 

R EAu ( 5 v 1 0 30 1 RIIJ 1,KU(J) ,Bfc TI 1JI ,BETU(J) ,$PLAIU 
WRITE! t, 1040 I Ri !J J ,RU< J) ,BbTI ( J) ,BE TlM JJ ,5PLNU 


MSP 1(3 1= SPLNJ 
NSP = N5PHJ) 
WRITE! 6, 1190 I 
READ I 5, 1030) 
WRITE! 6, 10401 
WRITE* 6, 1200) 
READ ! 5,10301 
10 WRITE! 6, 1040 ) 
WRITE! 6, 12101 


J 

(MS Pi l , J I ,1 =1 rNSPJ 
CM SP II « J ) 1 1 = I ,N3PI 
4 

i THSPd ,J) ,1 =1 ,NSPi 
i I H $P < l ,J) ,1 =i ,NSP1 


READ ( 5, 1030 ) 
WRITE! 6, 1040 1 
WRITE I 6, 1220) 
READ ( 3, 1030 ) 
WRITE! 6, 1040 1 
WRITE! 6, 1230 l 
READ I 5, 1030 I 
WRITE! 6, 1040 1 
WR ITE! 6, 12401 
READ I 5, 1 0 10 ) 
WRITE! E, 10201 


I MR I I ) i 
imi i i, 

( ft M SP (I 
<RMSP! I 

(BE SP U 
! B b SP ( 1 


-i ,NRSP1 
= i jNRSPJ 

,1 =1 ,NR$P) 
,1 =i ,NRSPI 


,1-1 
,i -1 


NKSP) 

NRSP) 


BtOAr,AANOR,ER SuR , STRF N , SLCRO, I MTV L# S OftV L 
BL0AT,AAN0R,EK$0R , STRF N , SLCftU , I NTVL,SURV L 
IF iHM-LE* 100.ANU.MB3! ,LB. 50 .A NO .Nft SP* LE. 50. AND. NS PH 1 l.LE. 52 
l.ANO.NSPI 1 1 ) *LE. bG.ANJ *NSPI (31 *LE. 5 0. AND. NSP 1 !4) ,LE.50 ) GO TO 


20 


WRITE 16, 12 501 
STOP 


CALCULATE MV ARRAY 


20 mi = CHORD! 1) /FLOAT! MbU-Mbl J 

IF! Mtki .OT .m -l2.AN0.rtb! .NE.M312) HMl = MLE (3 1 / FLOAT ! MB l 2-mi I 
mz = i. e30 

IF! MB i2.hlE.Md3 ) HH 2 = I HLc (31 -CH ORO II) ) / FLOAT ( MB 12 *MBU ) 

HM3 = CH3RUI 3)/FLUAT!MBU2-M&I2 ) 

IRMbQ .GT *m 1 2 * A N D # MB U * Nt . H6 02 ) HM3 = (CHORD ! 3 1 F MLE f 3 ) -CHURO! I ) )/ 
IFLJATi MBU2-MBU I 
MBUT = MINGI I 21 
C0MB3T = AM 1 N 1 (CHURO f 1 ) , MCE ! 3) ) 

00 30 IM^1,MB0T 
30 MVIIMI = FLOAT! IM-MBI l^HMi 
MBIT * MA K 0 i H60 ,MB 121 
COMBIT - AM A X 1 ! C HO RO ( 1 1 , M LE I 3 ) ) 

OD 40 IM=MBU T , MB 1 T 

40 MVIIM) = COMBOT+FLUATUM-HBUT) *HM2 
00 50 IM * MB I TjMM 

5 0 mv 4 in i = combi t+flOat i m-m i ti*h m3 
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C AL CUL AT t MISCELLANEOUS CONSTANTS 
N£R( U*Q 

ner( 2 mo 

PITCH * 2.*3*1415927/FLQA TtNBLl 
HI = P ITCH/fLUATiNB^n 
DTi R = HT/1UQO* 

miii = AHINli HMUHM2THM3J /lOUO. 

QVi II - 0 * 
fcJVt 2 ) = 1, 

6VI 3 I = -NTFLSP/wTFL 

avc = i.+a vt 3i 
naiMi- m i - 1 
MaiPi* m&i+i 

HBQM l * MbU~l 

MBUP 1= HtiU+i 

MBI2M 1~ MB 12*1 

MB12P1* MB12H 

HBQ2Ml~ M By 2*1 

MBU 2P 1 = M602+1 

WMM i = MH*1 

CP * AR/l GAM-1* )*GAM 

EXPQN- i./(GAM~ 1* * 

T W W “ 2**UMEGA/KTFL 

CP T IP - 2, ♦CP* TIP 

T GROG = 2 # * GAM* Aft / ( GA H + 1 * I 

CALL SPLihinHrt ,RMSP r NftSP,MVfMMtRMrSALI 
CALL SPL lWT[fKiBESP t NR$P * M V t MM P B E ,DBOMI 

CALCULATE GfcUMtTKlCAL CuNSTANTS 

CHDROI 2 i * CKlkDl II 
CHURUf 4) ” CHJKD( 31 
STGRC 2 ) * SlGRi 1 I 
STGRI4I * STGft(3I 
MLEI LI * 0* 

ML tU * 0, 

ML £t 4 i = MLE13J 
THLEi i I * 0* 

T bkL fcl 2 J = Pi TCH 

THLE1 41 - P I TCH* THLE I 3 ) 

KHi( II = RHi MBi I 
RM M 2) ^ RMI MB I * 

KM It 3 J = KMt mill 
RMlt 41 =* KMl MB 1 2 J 
RMUi n = KMiMbUI 
RMOi 2 I = RM(MBUl 
RMUI 3 1 = RMt MBU2I 
KMJ(4| - RHimOZi 

INITIAL I Z £ AKRATS 

m 60 i“ 1 1 2000 
U( i ) * U 
KU) ^ Op 
60 RB04 I I = Rri3 IP 
DO 70 IM“ 1 * 1 00 
DO 70 SUKF*lr4 
RHUHtH 1M, SURFI = RHOIP 
RHUVBt IM# SURF1 * RHOiP 
7 0 iTVt IM * SORE I = - 10000 
RETURN 

100 U FORMAT i IHH 
L010 FORMAT (1615) 

1020 FORMAT ilX t 16l7* 

1030 FORMAT (BFiQ.Si 
1040 FORMAT t IX , BG1 6* 7 1 
1100 FORMAT (BOH 
1 

1110 FURMAT ( 7X # 3HGAM,14X,2HAR,13X,3HTI Ptl2X,SHKH0iPtl2X,4HWTFLaiX*6H^ 
IT FL SP , 10X r SHOMEGA 1 12 Xt 3H0RP I 

1120 FORMAT t 6 Xi 5 Hb£TAl , 10 *, 5 H 6 E TAU , 11 X , 6 HCH 0 RQF , 1 1 X » 5 HST GRF ♦ ID X , 
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16HCBC1RDK* 10X,5HSTGRft * 12X * 4HMLE R , U X,5HTHLER) 

1 13 0 FORMAT I 41 H MB I W30 Mi i 2 MB 02 MM NBBi N8L NftSPJ 

1140 FORMAT 153Hl BLAUfc SURrACt 1 — UPPER SURFACE * FROM! BLADE) 

1150 FORMAT 4 5 3HL BLAUt SURFACE 2 — LOWER SURFACE - F HUNT BLADE} 

1160 format i52hl blade surface 3 — upper surface - rear blade} 

1170 FORMAT { 52hL BLADt SURFACE 4 — LOWER SURFACE - REAR BLADE) 

1180 FORMAT i 7X , 2 HR 1 , 1 1 » 12X , 2HRU , 1 1 *l2X,4HBETi ,11 , 11 X ,4HBET 0, 11 , 1 IX , 5H5 
1PLNG, i 1) 

1190 FuKMAT i ?X* 3HM SP * i 1 , 2X, 5HARRA Y) 

1200 FORMAT ( 7X,4HTH$P, tl, 2X*5H ARRAY) 

1210 FORMAT I 16HL MR ARRAY) 

1220 FORMAT (7X,11HKMSP ARRAY) 

1230 FORMAT <7X,llHBESP ARRAY) 

1240 FORMAT 152HL BLOAT A4NDK ERSOR STRFN SLCftD I NT V L SJR7U 
1250 FORMAT (4iHl MM ,N&B I ,NKSP * UR SOME SPLNU IS TOO LARGE) 

END 


SUBROUTINE P REGAL 

PRECAL CALCULATES ALL REQUIRED FIXED CONSTANTS 
COMMON SRw, i TfcK , 1 END ,LtR 121 ,N£R(2i 

COMMON /AUjUNO/ A 12000,4) , 0( 20001 ,K(200U I r RMOI230J ) 

COMMON /INP/ GAM, Aft ,TIP,KHOiP , WTFL * WTF LS P * UMEG A, ORF, 8£ T A I, BET AD, 
1H B I * M BL), 4 G 12 ,MB3 2 *MM * NdB I , NB L , NR 5P , MR 1 50 ) ,KMSP15Q) ,BESP(53 } ■ 
2BLUAT, AANUR, EKSJRtSTRFN, SLLftO , I H ? VL , SUK VL 
COMMON /CALCDjW Mb I M 1 , MB i P 1 . MB UMi , MB DPI , MB 1 2 Ml , MB! 2 PI * MB 02 Ml , 
IMBU 2P IfMrtrtl, HM1 *HM2*MM3 ,hT, 0TLR,DMLR , PITCH ,t P * EX PUN , T WW * CPT IP, 
2FGRJG, IBI, TBO , L AMbU A * T WL , I THIN, I TMA X , NT P , l MS (4 ) , BV ( 4) ,MV1103), 

31 VI iOl), 1 TV( 10C,4) *Tvl 100*4) ,0 TDMVUO0 ,4) , BET A V 1 13 0 ,4 ) * 

4MH1 10 0,4) ,DTUMH( 100,4) ,3fc TAM (1 00 ,4) »KMH 1100,41 *BEH{ 133 ,4), 

5RM1 100 )*B£< 100)*0aUMi 1U0) *SAU100) *AAA UCO) 

INTEGER BLOAT , A AND K, Eft SOK , STRF N , SLCRD, SURVL, A A TEMP ,S JRF *S JRFBV , 
IFIRST * UPPER* UPPRBVtSi ,sr,$Rw 
REAL K , X A X t L AM B LJ A , LM A X , HH , ML E * MR *MSL * MSP * M¥ , MV 1 Ml 
EXTERNAL BL 1 ,b L2, BL 3, i L4 

CALCULATE LAMBDA AND VI 

BETA 1 = BETA 1757.295779 
BETAO = BETAD/57. 295779 
T B I = SlNlaE FA M/CUSIBETAI) 

T BO * SlNlBt TAul/COSicieiAU) 

10 RHliT * RHO IP 

RBQVl * *TFl /BtlMB H/ft I TCH/COSIBfcTAl } /RFHMBl ) 

20 VI = RHJVI/RHJT 

LAMBDA = ft M( MB I J * 4 VI *5 1 N 1 BE TAI ) *■ QME GA*RM 1 MB! ) ) 

TT IP = l*-( V I** 2*2 **OMl 3A*LA MBDA— 1 UME GA*HM 1 MB1 ) f^2}/CPUP 
IF( TT IP.LE.O.) GO TU 30 
ft HUM B 1 » RHU IP* TT IP** EXPUN 

iFI ABSCHHUMB nRHDF) /RMUI P.LT*.O0O0Ol ) GO TO 40 
RHDT * RHUMB 1 
GU TO 20 

30 WTFL ~ WTFL/2. 

NER ( 2 1= NEH1 2> + l 
WR 1TE( 6* 1020 ) wTFl 
IF4NER4 2) .EU . 10) STOP 
GO TU 10 

40 VI = RRGVI/HHOMBI 

LAMBDA * kMIMBU*! V1*S1NI BE TAI ) + UMfcGA*RMl MSI) ) 

CALCULATE MAXIMUM VALUES FOR KHU*W AT LEADING AND TRAILING EDGE 

TWL = 2.*UMEuA*LAMBUA 

AA = 1 TWL-IJMfcGA*ftMiM3II J**2)/GPTIP 

TPP = T IP * ( 1 *-AA ) 

BB = F GRO G*TPP 
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TTIP * L.-Bt* /CPUP-AA 
KNOT a RHUIP*TTIP**EXPQN 
RHOtfHI = R HO T* SQR T f tJB I 
AA = I TWL— 4 UM£6A*RMlMi502) 1**2) /C PT I P 
TPP * T IPM 1 *-AA } 

6B = TGKOG*TPP 
TT IP * I.-Bd/CPTIP-AA 
KNOT = RMU 1P*TT1P**EXP0N 
ftHDHMO = RHJ T* SQR T I BB I 

CALCULATE VU AN D W~CkI?lCAL AT BLADE LEADING AND TRAILING EDGE 

KHUVO - WTFL /BEIMB02J /PUTH/COStBETAGl /RHIMBQ2 J 
WHOM B 2 = KHJIP 

TWLMR = T ML- ( JHE GA+RH I MBQ2 I 1**2 
LERI L I =1 

OEN5TY CALL NO * I 

CALL DENS T Yt KHU VU , KH J MB2 ,VU,TWLMR, CP TI P,£XPGN,RHGIP, GAM, AR, T IP1 
HCRI = $UKT< TGRDG*TiP*U,~t TWL-lOMLGA*RMItt6m **2 J /CRT 1PU 
RCKO = SgRTI T&kOG* U P * (1 - - I T HL- I UHt GA*RHI M802 I I **2 J/CPT IP) I 

CALCULATE BETA CORRECTED TO BOUNDARY A-N AND G-H 

T mi. MR = TWL- (UH£GA*KM[ 1H**2 
RHOl = RHUMB I 
TBII = 1.E20 

50 TBIT ^ ITBi/BEtHSl l*R HOI /RHUMB I + GMEb A* I RM < MBI) **2-RMl 1 >**2»*KHG1 
L/rfTH_*PITCH)*Bfct 1 1 
£ FI ABSCTB I i- TBI T }.LI, ,000011 G O TO t*Q 
TBIL = TBIT 

KHUVl - riTH. /P ITLH^SOR U l. + TBI 1**2 I /BE (i l/RMI U 
L ERI I I =2 

OENSTY CALL NO* 2 

CALL DENS I Y i R HO VI ,RHU1. AA ,T HLMR t G PTI P ,E XPUN , KHU I P r GAM, AR, T 1P| 

GO TO 50 
60 iBi = TBIT 

BTAIN * A TAN f T8J 1*57,295779 
TwLHR = TWL- 4UH£GA*RMI MMH**2 
RHOMM = RK0HB2 
T0QM = 1.E20 

7 t> T 6Q T * ( TBD/&E4MaG2l*KHUtfM/KHUMB2+UNfcG A*I RMI MB02 J**2-RM( MM )*#2 1 + 
LRHQMM/*TFL*P ITCH I *0E 4 HM1 
IF I AB 51 T BUM -TBO T j *L T » *00001 1 GO TO BO 
T&UM - I BUT 

RHDVD * *TFL /F I TCH*SURTIi ■ + TB0M**2l /BE IMMI /RH(MM) 

LERI 1) =3 

QENSTY CALL NO, 3 

CALL DENS I Y iRHUYUiRrlUrtM, AA, Tv4LMk,CPTl P,EXPON T KHUI P r GAM, AR, T IP 1 
GO TO 70 
BO T BO = IBJT 

BT A OUT = AlANl T60l*57. 2957/9 

CALCULATE TV, ITV, IV, UTDMV, AND BET AV ARRAYS 

ITNIN = □ 

ITMAX = NBBI-I 
C TV, ITV, AND DTOrtV UN BLADES 
DO 90 IH-rtBIt«BU 
LERI 2 1 = 1 

C BL CD CALL NO. I 

CALL BL II M VI IM I, TVI IH, 11 ,0 TD MV il H, U , I NM 

IT V I IK , L J -= I NT U TVI I M , 1 1 *0 TLRi /fill 

If l TVI IM, II .GT.-OTLR I l TVI i M,U =1 TVIlMtil+l 

ITM IN* HINDI 1 TM f N , I TV { IM « 1 1 1 

LERI 21 =2 

c blcd Call nj . 2 

CALL BL2I MV4 4M l t I VUM,2) ,0 TD MV 4 1 M,2 ) , I NF i 
IT VI IM, 21= 1NTH TVI 1M,2>-DTLR1 /HTl 
IF I TVI IM , 2) .LI.DTLR I I T V I I M ,Z I = IT VI I M ,2 I -l 
90 ITMAX^ HAXOI ITMAX, I TV< 

DO llO IM=MB 12,M&u2 
LERI 21-3 
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C BLCU CALL NJ , 3 

CALL OLBIHVi iHJtTVCIM, 31 ,JTUMVU M,3l ,1 NF) 

ITV< INTIUVUN, 3 1 + 0 TLR1 /H Tl 

if ([Vi^ f 3KGn-urui irvuM ( 3>=imi«.3m 
IF llrt.bT.MdUl CO TO ICO 
TV* IH,,3) = TV* Irt,3I+Fl TCH 
100 IT rt I N = NlNOUTHtN , 1TV{ IM»3) 1 
L£R( 2)^ 

C BLLD L AlL NJ, 4 

CALL 8L4* HV< IMIrTVUM ,41 ,QTOMVl IH v 4t , I NF I 
1TVI IM,4l = lNTUTVUM f 4}-QTLRl/HT) 

IF f TV( W r 4J .L T.UTLKl 1 T V * IH ,4 1 = I T V ( I H ,4 ) -1 
L10 I TMAX = HAXO* I TNAX , I T V 1 i H *4 M 
C I TV AND IV UPSTREAH OF FRONT BLADE 
FIRST ^ 0 
LAST = IHBB I- 1 
DO 120 IM = 1* HB I H 1 
IV* 1M > = NBB 1*( IH-ll + l 
IT V ( 1M , 1 J = FIRST 
2 20 IT VC IM,2J = LAST 
C I TV &ETWfctN FRONT AND REAR BLADES 
IF <H6UP1 -GT -M8| 2H11 \,0 TO L40 
LAST* ITVtWB 12,4 J 
FIRST - LAST+l-NBGI 
00 130 m=MB0PI,MBI2Ml 
1 T V ( 1H,31~ FIRST 
130 ITV t IM , 4 1 * LAST 

ITHiN = H INO( 1THIN,I TVIM3UP1 ,3H 
C ITV DOWN STREAM oF REAR BLADE 
140 L Ai T - ITV(MBU2 r 4) 

FIRST = LAST+l-MbBI 

00 150 IH=MdU2Pl,KM 
iT VI IM,3l = F IK ST 

150 ITV{ IH , 4 1 = LAST 

liniH = M IHOf ITHiN r lTVIMM,3M 
C FINISH CALCULATING IV ARRAV 
IViHdM * NBOI^M&lMUl 
MBOT * H 1NQIHBU ,HB I 2M IT 
IFIMtH .GT .MdUTi GO TO 165 
DO 160 lN=M8i f HBUT 

160 IVC1H + U = lVUMlMTV<lrt v £|-iTVUM t Ll + L 
165 IF* H8 1 2 *GT *M bO 1 GU TO 180 
DO 1 7 0 IM =Md I 2 , H50 

170 IV* IM + 1T - 1 Vi IMj + I TV* IH,2J-1 T VI IM r 3l* I TV 1 I 1 -ITVUM, 1T+2-NHBI 
180 DO 190 IM=M8QP 1 ,MM 

190 IV*IH + 1T * IVI IM J + ITVI I«,4)“lTVl 1H,3) + 1 
C BETAV AKRAV 

DO 200 SUft F = 1 , 2 
DJ 200 IM=HBI,H83 

200 8ETAV* IM, SUKFi * A TAeW D I J HV U M , SUkF 1 *RMU H)i *5 7 .29 5 7 79 
DU 210 SURF = 3, 4 
DO 210 I2,HBU2 

210 BET A V I IH , SUR F I = A TA N * 0 OH V i 1 M , SURF 1 *RM [ l Mil *5 1 * 29 5 7 79 
NIP * iVtHHl +NBBI-1 

WR I T b 1 Erl 030 1 VI .AHUwfll , *ICiU ,B TA In , VO , RHU WHO , WCRO, BT AUDI 
MR [ TE { t, 10401 P I TC H , H T t H H 1 ,H M2 ,H M3 
m 1 T El 6, 1050 I I THIN, I TMAX, LAMBDA ,NI P 
WR I TE ( 6, 10601 ( SURF ,8 V ( SURF ) ,SURF = 1 ,4* 

IF * 6L OAT »L t » Q ) Gu TO 230 

FIRST * MB I 

LAST * MGU 

hR lit (6, 10701 

DU 220 SURF= 1,3,2 

1 = SURF+1 

WRITE (6tl08Gl $URF,i , (MVUHT , TVII M, SURF I ,DTDMV( IM,SJRFT * 

ITV* IM, 1 1, DTD MV ( IH, 11 1 IM=FlRST,LASn 
FIRST = HBI2 
220 LAST = mo 2 

WRITE 16, 10901 t IM'HVUMJ ,RHU HI rSAL(IM) f 6E(l Ml ,D80H( I Hi , J M= 1 » MM 1 
230 CONTINUE 
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C 

€ CALCULATE MH A NO DTOMH ARRAYS 
C 

ITO = IT V l L* li 
MftfS - 1 

msi i I = I 

MHt 1, 1 J = 0* 

UTDMHt I* 1 I = 1 »E 1 0 
LEKi 2 MS 

C dLCO AND ROUT I VIA HHUKiZ) CALL NO* 3 

CALL MhORl/IHVrlTVatn , MB Q t l T QpBT *DT LR ti> p I HS (U » MHi i, II, 

iPrOMHC It 1 I f HR T S I 

IF ( I T VI M8Up I ) - f TVfMBU t Z ) f N3BI * ME. 2} GO TO 240 
W$L - IMS! 1 M-L 
MHi IMSL* i I = HV(HBU) 

DTQHNC i,4 SL t 1 i = ~I*£iO 
m$i U = 1MSL 
240 1HSI2) = 0 
MRT$ = 1 
L ER{ 2 f =6 

C BL CD A HD ROOT (VIA MHURIZ) CALL NO* 6 

CALL MEOR U(MVtlTV!l v 2l ,BL2,MBI ,MQ f l ra,HT,OTLft,l t I MS ( 2 I * MHI 1, 2 I, 
IOJD Mhl 1* 2JtHRTS| 

IHSI3I = 0 

If 4 I T V(MB 12 ,i HI TVf MB !2,4) *NE*2 .AND. 

IITVIM6 12* 4H iTVIM8l2p3i.NE.NBBI-2l GO TO Zb 0 
MKTS * 1 
IrtS 4 3 I - l 
HHi 1, 3 I =* MViM8I2) 

DTDrtHt 1* 3) - I.EiO 
250 lEKI 2) -7 

C atCO AND ROOT (VIA MHJRiZ) CALL NU. 7 

CALL MHJK IZ( Mv t I TVU ,3) ,OL3*H6i2 trtGO*I TO, NT ,DT LR ,3 * t MS ( 3 ) P HHi i , 3 } , 
1QT DM Hi l*3|,riKTS) 

H BUT -* MAXO i MaO # HB 1 2 1 
LERI 2MB 

C BLCU AND ROOT (VIA MHUfUZl CALL NO. 8 

CALL rt K)R IZiHVtlTVilpBI , BL3 . MOOT *M602 * 1 TO ,HT ( UT L R * D , I MS l 3 J * 

LHH( i, 3 It OTD* MI 1,31 t MRTSI 
IF ( I T V( MBG2, 3M i TVIHBU2 ,41 + N6 6I *NE* 21 GO TO 260 
IHSL = IMSI3HI 
HHi M$L, 31 * MVf MBU2I 
UTOHH( 1HSL*3I = -1.EI0 
IM Si 3 1 = I MSL 
260 IHSI4J = 0 

if ( IT VI MB 12 til- ITVIMBI 2»4KE0.2.0R. I T V i MB 1 2 , 4 ) - l T V ( HB 12 , 3 1 . EQ * 
1NBBI-2J HR TS = 1 
LEKI 2M9 

C BL CD AND RJJ I I VIA MHURIZ) CALL NO* 9 

CALL MKiR IZi HVplTViltAI , SL4 , MB I 2 * HBG2 1 1 TO ,HT ,01 LR * 1 *IHS(4), 

£HHf l f 4 JtDTDHHi it 41 * HR T SI 
i * HAXOi iHSIlMIMSiZl ,IHS(3I tlHSi4M 
IF i ULt .LOO I GU TO 2 70 
rfRITULtiiOO* t 
STOP 

CORRECT I TV ARkAY FOR Q VtRLA PPI NS PORTION OF BLADE SURFACE 3 

270 IF (H0 I2.GI.MBCl I GO TO 2 90 
00 280 IH -MB 12* HB0 
260 ilVC I M f J } * nvilMpSMNSBI 
290 iFiBLOAr.LE.OI GU TO 300 

WRITE (6* 11 1C) t 1M»I ViiMt ,if TVUM,$URF| .SURF^l *4) v !ri»l »MM) 

CALCULATE Rrth* 8EB* AND fl E TAH ARRAYS 

300 IFitiLOAT.GT.OI WRITE16*U20) 

OO 320 SURF - 1*4 

CALL SPL lNTtMK,RHSP ( NRSP,MHU,SURF) .TMSISURF) , RMH U , SO RF I , A AA ) 

CALL SPLINT! HR t BE SP , NR SP * MH U p SURF I t IHSiSURF» , BE H i I P S J RF i • AAA i 
IMSS * IMS! SURF) 

If( IMSS.LT.1 ) GO TO 32 0 
00 310 IKS = It 1 M 5 5 
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310 EJETAHi IMS 9 SURF } = A TA N ( J TO MH U r! S , SURF J *KMH (IH S ,SURF H *57 .2 95 7 7 9 
if IBlCAT.GT.D* tfht | TE 16,11307 SU RF # I MH ( I M , SU RF ) , RHB 1 1 M f $ J R F > t 
1BEHI (M,SURF* .dtMHUM, SURF J , D TOMH I i H r SURE } f l M= 1 t IMSS) 

320 CONTINUE 

IF t BL CAT ,Lfc ,01 GU TO 3*0 
WRITE I6 f 1140* 

IT = ITMIN 

330 IF I I T ,GT . ITMAXi Gu TO 340 
TH = FLOAT! i IH'HT 
kkHL 16,10107 1T,TH 
IT * It+i 
GO TO 3 30 

340 IHN£F 4LE.200O} GO TU 350 
m ITE1 t , 1 1 50 ) 

STOP 

350 rtR iTE i 6 r 100 0* 

RETURN 

1000 FORMAT I INI i 

1010 FORMAT 1 4X ( I 4, G 16*51 

1020 FORMATtGOHLlNPUT WEIGHT FLOW IwTFU IS TOO LARGE AT BLADE LEADING 
1EUGE/16H WTFL REDUCED TU,Gl4.6| 

1030 FORMAT C 1 H 1/ // / / 24X, 10HFKE L S IRE A M, B X T l 3HHAXI Mti M VALUE, 

1 7X 1 SriCR I T ICAL, 30X, 14H&ETA C uKREC TE D/25 X ,SH VE LOG I T Y , 1 0 X ,9HF OR RHU^rf 
l, 10X, 3FVEL0L I TV, 3iX, ilHTO E UUND A RV / 1 X , I 7HLE AO i NG EDGE B-M,3G1B,5, 
312X* IZhBOUNDAHY A-N ,GlfcU5/iX,l 7H TRAI LING EDGE F- 1 , 3G 13 * 5 r 1 2X * 
412HBOUNOARY G- H , G 18« 57 

1040 FORMAT I / / // / 5X ( 2SHC A LC UL A TtO PROGRAM GUNS TANTS//3 X,5HPITXH, 13X, 
12HHT, 13X, 3HHM1, 1 3X , 3HHM2 , 1 3 X , 3HH M3 / 1 X ,5G i 5 * 7 7 

1050 FORMAT U 5X, 5H I TM IN , 10 X , SHI TMA X/ *X * I 5 1 1 0 X , [ 5/ / 5 X *5 Hi AM QUA/ IX , G16*7 
1/// 5X t 33HNUMBER OF INTERIOR MESH POINTS = ,15* 

1060 FORMAT l tftffi 5X,23HSURFAGt BOUNDARY VALUES/ / 5 X , 7MSU RF AC E, 7X #2 H8V 
1/1 3X, I 4, 4X,F 10 # 5 J 1 

1070 FORMAT I i H I, 6X * 62HdLADt DATA AT INTERSECTI CN5 OF VERTICAL MESH L IN 
1ES WITH BLADES* 

lOdD FORMAT I 1HL, 22K ,13H BLADE SURF ACE ,1 2 U 5X , 1 3HBL ADE SO Rf AC E * 1 2/ 7X , 

1 1HM ,l4X,2HTV,llX f 5MJTJMV,12 X*2HTV ,1 i X ,5 HDTO MV/ ( 5 G 1 5 • 5 I * 

1090 FORMAT l LH l , 13X , 44HSTRE A M SHEET COORDINATES AND TH 1 C RN&SS TABLE / 

1 2X, 2HIM, 7X , 1HM, 14X,IHR i!3X,3H SAL 1 1 3 X , L HB , 1 2 X ,5H OB/ DM/ { LX, 13, 

2 5G15.5J7 

1100 FORMAT I 34HLUNE OF THE MH ARRAYS IS TGU LARGE/ 7H IT HAS, 15, 8H POi 
1NTSI 

1110 FORMAT I4H1 IM,9X, 8HI V ARRA Yi32X,9HE TV A RR A Y/ 3 8 X ,5 H BL A uE/ 37X , 7 HS J R 
IF ACE, 3X 1 1 H L, 5X, IH2 , 5X , IH 3 , 5 X , l H4/3 9X , 3 HNU< / i L X , I 3 , 5 X , I 13 , Z5X , 

24{ 14* 2X7) * 

1120 FORMAT 16 7H1M COORDINATES OF INTERSECTIONS OF HORIZONTAL MESH LINE 
IS WITH BLADE! 

1130 FORMAT ( 2 5HL MH ARRAY - BLADE SURFACE , 1 2/ / i 5X ,2 HMH , 19X t 3HRMH , 19X , 

1 3HBEH, 1BX , 5HBETAH,17X t 5HD T3 MH / 1 5G22 . 4 1 ) 

1140 FORMAT W3H1THETA COORDINATES OF HORIZONTAL MESH L INES//6X , 2Ht T , 
15X, 5HTFEIA * 

1150 FORMA M48HL f HL NUMBEk OF INTERIOR MESH POINTS EXCEEDS 2030 I 
END 


non o fir rt n rr o n n r noon 


SUBROUTINE COE F 


COEF CALCULATES MNlTt D I FFfcKENC E COEFFICIENTS, A, AND CONSTANTS, K • 
AT AIL UNKNOWN MESH POINTS FOR THE ENTIRE REGION 

COMMON $KW, £TtR*ItNO,LER(2J ,NER(2) 

COMMON /A UK l HJ / A 1 2000 ,4* , Ut 2 000 1 * K (2000 I , RH Ui 23 03 f 
COMMON /INP/ GAM i AR , T I P ,RHUI P , W TFL , RTF LSP , OMEGA, QftF*8ETAI , BET AO, 
lMBI,M0U,MBl2,MBO2,Mrt*NBBI , NS L ,NRSP , MR ISO I ,RMSP(50I *&ESPt bO J, 
2BLDAT, AANDK, ERSDR, STRPN, SLCKD , 1 NT VL* SURVL 
COMMON /CALCU HJ MB I M 1 , MBi P L f MB GM1 , MBOP1 , tf BI2 Ml , MBI 2 Pi T MSG2M1 
1MBO 2P1,MMM1, HMl,HM2,HM3#HT # OTLR f OHLR ,Pl TCH *CP , EXPON ,T RW , CPT iE 
2TGRQG, IB|, rao, LAMBDA rTWLtlTMlN, I T MA X , N 1 P * I MS (41 , 8V I A ) , MV ( IDO > * 
3IVI 101 n t T V( LOO, A ) r T Vi 100,41 ,DTDHVUOO,4I * BET AV (130 ,4 ) * 

4HHt 100,4I,DTOMHUOO t 4l,SETAHllOO,4| ,RMH (130,41 ,BEH(103,4I, 

5RHI 10QI,BEI lGOJ,oaOMUOO> , SA L ( 1 001 ,AAAUOOi 
COMMON /HR BA AK/ H ( 4 I ,R (4 1 *B I4l * KAK ( 4T , KA (4 I , 1 H (4 J , RZ , BZ 
INTEGER bt DA T , A AMOK, ER SDR , STRF N , SLCRD , SURV t,AATE MP *SJ RF *S DRF BV , 

IF IR$T , UP PER, UPPRBVtSlf ST , SR W 
REAL K i KAK , LAMBDA , LMA K ,MHl , MLE ,MR ,MSL ,MSP *MV* MV I Ml 
INITIALIZE ARRAVS 
ITtK * ITtR+l 
iHt U - 1 
DO 10 1-2,4 

u> ini n = 0 

I FI I TV (MB 12, 3|“ItVi MB I 2 , 4 J *E 0* 2 J £ H t B I = 1 

IFf ITVIMd 12, 4I-iTV|Mai2,3i*EW- NBBI-2, AND.M8I2* N£*MBUP1 I IH(3I = 1 
I NCOM PRESS! BL E CASE 

1F( GAM *NE *1* b*UR * AR *Nt * 1 OuQ* • OR# T1 P* NE - 1 • E6) GO TO 2D 
I END st i 
GU TO 4Q 

ADJUSTMENT JF PRINTING L ON TR QL VARIABLES 
20 1 F I I T ER *N E # 1 «AN0 * I T LK - Nfc * 2 J GU TO 30 
AANUK * AANUK-l 
ERSOR * ERS3K-1 
STRFN = STKFN- 1 
SLCRU = SLCRlF- l 
IN I VL = INTVL-I 
SURVL = SURVL- 1 
30 IF{ ILN D.V E *0 > GO TO 40 
A AN UK = A A NO R + 2 
ERSUR = ERS3R+2 
STrFN * STKFM+2 
Si CRD - SLCRD+2 
SNTVL * INTVL+2 
SURVL * 5UKVL+2 

FIRST VERTICAL MESH LINE 

4Q UO SO IP s 1 ,N 88 I 
At IP, I J - 0- 
Ai IP, 2 I = 0* 

AT IP, 3 J * 0* 

At IP, 4 J ^ 1* 

SOKtlPJ * HMl*r01/P!TCH/RM4iJ 

yPSlREAM UF BLADES* EXCEPT FUR FIRST VERTICAL MESH LINE 

lFt2-GT,MBIMII GU TO 70 
DO 60 IM - 2, MB I M 1 
60 CALL COtFPI I M , 1,21 

BETWEEN FRONT BLADES 

TO MBOT ^ MIN0(MB0,MBI2MU 
IFLMBI *GT ,MBOT| GU TU 90 
DO BO IM -MB I ,MBO T 
80 CALL CUEFBbi lM,i*2*n 
9U !FtM&I2*GT.M8U» GO TO 110 

OVERLAPPING CASE 

00 100 1H -MB 1 2 , MESG 
CALL CUEFBBi IM , 1,4, U 


x * 
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Nil III III III III! III III^M III llllllllll 1 1 


16 0 CALL COEFBBi IM, 3, 2, 41 
60 TO lit) 

NON - OVERLAPPING CASE 

110 1FIM0UP1.CT.MBI2.111 60 TO 130 
00 120 IM =M6QP 1 1 MB 1 2M1 

120 CALL COEFPUM, 3,41 

B£T»fctN REAR BLADES 

110 MBIT » MAXOt MB I 2.MBUP 1 1 

IMMBU.GT.M3U2I 60 TO 130 
00 1 AO IM=MB IT , MBD2 

140 CALL COEFBBI IM,3,4,3l 

DOWNSTREAM OF BLADES EXCEPT FOR FINAL MESH LINE 

ISO IFIMB0 2P1.6T.MMM 11 GU TO 170 
OO 160 IM =MB02P l.MMMl 

160 CALL COEFP < 1 M, 3, 41 

FINAL VERTICAL MESH LINE 

170 I VMM = IVIMM1 

00 1B0 IP -iVMM t NIP 
Al IP. 11 = 0. 

A< IP, 21 * 0. 

A( IP, 31 « 1. 

Al IP, 4 1 * 0. 

180 Kt IP} = 'HM3*TB0/P I TCH/RM1MM1 

TAX 6 CARE UP POINTS ADJACENT TO Bt AND CAStS WHEN POINTS J,C, E, UR F 

ARE GRID POINTS 

POINT B 

IP = IVlMBHll 
At IP, 41 = 0. 

POINT J 

IFt ITVIM8I2, 31-lTV(MBI2,41.Nt.21 GO TO 190 

IT = I IVt MB 1 2, 4 1 F 1 

IP « IPF( MB 1 2M I , i T 1 

XI IP 1 ■ R( IP }*A< IP , 41*8 Vl4l 

A< IP, 4 1 * 0, 

POINT C 

190 IFt ITVIMB0,11-1TV(M80,21*NBB1 .NE.2I GO TO 200 
IT = I TVt MBJ , 1 1“ 1 
IP = IPFtMBOP 1,IT1 
At IP, 31 = 0. 

POINT E 

200 IFt ITVtMB 12,4)-ITVtMB12,31.Nt.N8BI -2. Oft. MB12. EO. M0QP1 1 GOTO 213 
IP = 1V(9BI2MH 
K( IP I * Kt IP i+Al IP,4I*BV131 
At IP, 41 = 0. 

POINT F 

210 IF (t ITVtMe02,31-ITV(MBO2,4lFNBBl.NE.2l.ANO. UTVIMBU2,3I 
1-ITV1MB02.41 .Nfc.211 GO TO 220 
IP * 1 Vt M BO 2P 1 1 
XtiPI * K t IP l*At IP »31*BVt 3 1 
At IP, 3 1 * 0. 

LINE X-L AND LINE TO RIGHT OF C-0 

220 FIRST = MAXOt ITVIMBOPI .31FNBBI ,1 TVtMBO.31 I 
IPKL » IPFtHBU.F 1KST1 
IPL = IPKL*-| TV<MB0,21-FlftSI 
IPCO = IPFtMBOP l.FIKST-NBBl) 

230 IFt IPXL.6T.IPL 1 KEIORN 

Kt IPXL 1 = Kt iPKLt*AUPKl,4l 

Xt IPC01 = Xt IPC01-AtlPCD,3l 

IPXL = IPKL + l 

IPCO = IPCO* l 

GO TO 230 

END 
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SUBROUTINE C ULFBB I IM, UPPE R , L UWE R »UPPRB VI 

CUEFBb CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K 
ALONG ALL VERTICAL ME5 H LINES WHICH INTERSECT BLADES 


COHMUH /A UKR HO / A 12OO0 ,4K U12Q00I , K {20QOJ , RH Gl 20 jj 1 
COMMON /InP/ GAM,AR ,TI P i RHUI P , W TFL , WTF L$P , GMEG A , ORE , BET A I , BET AO, 
IN 01 ^IBU, MB 12 ,1*002^, N0Ul , N0L , NR SR , HR i 50 1 , RMS P (501 , BE5 PESO I , 
2BLDAT, AANDR, Ek SOR , STRF N , SLCkD , I N TVL , SURVL 
CUMHUn /CALC UN/ RBlHWMtllPl , MBUM! , MB QP1 , KB 12 Hi , MBi 2 PI , rt0O2Ml , 
1MB0 2P i,MMMl,HMKrtM2,HM3,riT,UTLrt ,OMLR,PITCH ,C P , EX PQN ,T W W , CPT IP, 
2TGRUG, ft* I t TdUf LAMBDA , THL , 1 TNI N , i TMA X ,NI P ,1 MS ( A * ,aV(4» ,MVU0O 1, 

31 Vt 101 U ITVI 100,41 ,TVllO 0,41 ,0TDMV(10O,4» , BET A V (UP , 4 > , 

4HH< 10 0,4KUTDMH1 1U0,4) ,BETAB(1 GO ,41 »RMH(1Q0.4I ,8EHU33 ,41, 

5RM( iUO ),BEI lOOi,D0UH(lOOI ,SALUOO> ,AAA(100I 
COMMON /HKBAAK/ H(4i , R ( 4 J , B I 4 > , KA K ( 4 ) , KA { 4 I ,1HI4),RZ T BZ 
INTEGER BLOAT, AANOK,bRSUR, STRF N , SLCRD , SUR V L , A ATEMP ,SU RF , 5U RF BV , 
IF IKSI , UPPER, UPPKBV, Si, ST, SRW 

REAL K ,RAK,L AMBDA,LMAX,tfH,MLE ,MR ,HSL , MSP , M V, MV I Ml 
I F I IT V I IM , UP Pfck KG T. I TVUM,IUW£ KM RETURN 
IT VU = ITV(IM,UPPEkI 
IT VL -» ITVl IM,LuWER> 

II “ ITVU - l 
IPU IPF( IM, ITVU I 
IPL = IPU* IT VL- ITVU 
DU 90 IP = IPU, IPL 
IT = Il+L 

CALL HKBt IN, IT, IP > 

00 10 1= 1,4 

RAM 11 = 0, 

10 KA< II =0 

FIX HR B VALUES FUR LINES C-D AND K-L 
1F( 1M.NE,M0JPU GU TO 20 
{FI IT «GEiITv(1M-1,1M GU TO 20 
HP 3 = IPFI IM-K IT+NBftl I 
R1 31 = RHUI i P 5 I 
20 IF! IM + 1.NE ,NBUP1 ) GO TU oO 

1 Pi MB 12P1-MGUP 1) 30,30,4 0 
30 IFI IT- IT V I EM , 3 J I 60, SO, 50 

40 1 Ft IT *LE * ITVI iM+l,4| I GO 10 60 
50 IP 4 = IPF( IM + 1,1 T-NBBII 
Rl 4 I = RHUl I P4» 

FIX HRB VALUES FUR CASES WHERE MESH LINES INTERSECT BLADES 
60 IF UT -ED* IT V I 1 M , UP PER 1 I L A L L B} M Y12 tl , 1 M, I T ,U PP ER, J PPRBV > 

IF III *EU # ITVUM,L3Wtkn CALL BD RYl 2 12 , 1 M, I T , L QW E R, L Oa EfU 
ITVM1 * i fvt irt-1, UPPER I 
ITVPi = ITVI IM+1, UPPER J 

IF I IM ,tO *MBU, AND. UPPER, EU.3I 1 T VP I = 1TVPH-NBBI 
IF i IM .Ej.MdLJP 1 I ITVMi - 1 TVMi-NBBl 
IF I IT ,L T * I T VM 1 J CALL BDR V34 (3 , I M, UPPER,UPPRBVi 
IF C IT -LT.iTVPlI CALL BOR Y34 ( 4 , 1 M, UP PE K , UPPR BV I 
1F1 1M .EU,MbI 2.AND.LUWLR.EU. 4) GO TU 10 

I F t I M * EQ « M GO P 1* A ND * L U WE ft * E U , 4- AND* MB 1 2 ♦ GK MS 01 GO TO 70 
IF ( t T .GT * IT V< IM-1 , LOWER I I CALL BO R Y34 i 3 , 1 M , LOWER, L OW fc R \ 

70 IH 1M ,£U*MBO *AND,LUWER*EU.2I GO TO HO 

IF I IT .GT . ITV( 1M4-1 ,LUHER1 ) CALL BDft V34 14 , 1 M, LOWE R, L OW E fU 
COMPUTE A AND A CUtPF IC I L N TS 
SO CALL A AK i IM, IP I 
DO 90 1=1,4 

Al IP I = K! IP l+KAKUDA t IP, I I 
9 0 iFtXAl IKEU.l) AUP r tl = 0* 

RETURN 


COEFP CALCULATES FINITE DIFFERENCE COEFFICIENTS, A, AND CONSTANTS, K, 
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ALONG ALL VfcRTICAL HE SH LINES WHICH OU NOT I NTERSECT 6LADES 

ENTRY CDEFPI IM, UPPER, LOWER) 

IFVU = J T V I IN, UPPER I 
IT VL - I T V I I H i LOWER J 
IT = IIVO-1 
IPO = 1V1 INI 
I PL = IV( IM+II-1 
OO 10 0 1P~IPU,IPL 
IT = IT+1 

CALL NRB1 IM, IT, IP) 

IF I I T *EQ ■ I T VO i Kill = RHOUPL) 

IF IIT.EU.HVL) R1Z) = RriUUPUl 
IH IM.Nfc,HBJP 1) GO TO 100 
IF! ir*GE*!TY(lM-l,lM GO To 100 
IP 3 = IP F I Irt-1, I T+NBO I ) 

K( 31 = Rrt j ( IP3) 

100 CALL AAM IN, IP ) 

M 1PL I = Kl I PL >+Al I P L ,21 
Kt IPO) = Ki IPUI-A i lPU,il 
RETURN 
END 


SUBROUTINE HRd< IN, IT, IP) 

HKB CALCULAELS He SH SPACING, H, DENSITIES, RZ AND R, AT GIVEN ANO 
ADJACENT POINTS, AND SIR t AH SHEET THICKNESSES, BZ AND B , AT GIVEN 
AND ADJACENT P3 ENTS 

COHHJN /AUKKHG/ A 1 2000 ,4 ) ■ Ul 2000 ) iK 12000 ) , RH 0 1 2QU3 } 

COMMON /LALCUVZ HB I H 1 , m l P l , MB GM1 , MBOPL , MB! l Ml , HBI l PI , M002MI , 
1MBU 2P l, MMM 1, HMl ( HM2,HH3tHT,L)TLR ,DMLR ,PI TCH ,C P, E K PON, T WW , CPT IP, 
2TGRJG, 131 , TdU, LAMBDA , T WL , 1 THIN, I TMAX ,NI P , I MS 1 4 I ,6V 14 I ,MVUOO ), 
31 Vi 101 ), I T VI 100,4) , TV( 100,4} ,0 TDMVUOO ,4) , BET A V 1130,41, 

4HHt 100,41 , Of DHH( 100,4} ,0 E TAH { I 00 ,4 } ,KMHU00,4) ,b£H flJD ,4 > , 

5RHI 100 1,0 ti 1 00 *,030*11 100) , SAL U 00) , A AA 11 OO } 

COMMON /HKdAAK/ H I 4 } , ft I 4) ,0 1 4 » , KAK 1 4 J , KA ( 4 > , I H (4 ) , RZ , BZ 
INTEGER BLOAT , A A N J K , e R SDR , S TRF N , SLC KD , Sllft V L , A A T b HP ,S J RF ,SU RF BY , 
IF IK ST , UPPER, OP PR B V , Si , ST , SR W 

REAL K,KAK,l AMBDA , LHA X ,MH , HLt , MR ,MSL ,M5P , M V , MV I Mi 
hi 1 )^ HT*RMI IM ) 

HI 2 1 - FT* RM ( IH) 
hi 3)- MV1 I« I-MVI IH-l) 

Hi 4) = M VI IM+ ll-MVl IM) 
iU = RhOI IP ) 

IP 3 = IP F I IH-l, IT) 

IP 4 = |PF( IH +1, I n 
R( 1 ) 35 RHO I IP- 1 ) 

RI 2 ) s RrtOi IP + 1 ) 

R< 31 * RHOl I P 3 ) 

RI 4 I = RHUi IP4J 
BZ = Bfcl IH ) 

BC 3*= Bel IN- 1) 

BI4J- BE! iMM) 

RETURN 

END 


noon 


SU&RUU T IN E AAAI IM, IP) 


A AK CALCULATES FINITE DIFFERENCE COEFFICIENTS, A » AND CONSTANT, K, 
AT A SINGLE MESH POINT 

COMMON /AUKK HO/ A (2000,4) ,Uf 20001 t K ( 20001 ,RriUC2D33 ) 

COMMON /C ACCOM/ M3 I M 1 , MB I P 1 , MB OMi , MBOPi , Mb 1 2 MI , MSI 2 PL , MB 02 Ml , 
IMB0 2P ItMMM i, HM L ,hH2,HM 3 , HT r Q TLR ,DM LR , PITCH ,C P, EX PON ,Taf tf , CPT I Pi 
2TGR0G, la I, la 0 , LAMBDA , THL,I THIN, I TMAX f NIP,I MS C A I , BV ( A } v MV l IDO I , 
31 VI 101 li ITV( 100,41 , TV ( LOO, 4) tOTOMVUOO ,41 , BET AV 1 13 0 ,41 , 

4MH< lUOrA) ,OU)MH( 100,41 ,3E TAH I LOO ,41 ,RMHtl00,4> , BEN CIO 3, 4), 

5km t 100 ),OEI 100) ,I>aOMl lUO) , SAL 1 100) , AAA UOOl 
COMMON /MKBAAK/ H i 4 ) , R i 4 ) , 0 C 4) , KAK 44 ) ,.KA < 4 ) , I H 14 I , RZ , B2 
INTEGER BLOAT,AANDK,ER 3UR, STKFN, SLCROrSURVLiAA TEMP, S JR F, SURF 0V , 
IF IK ST, bPPERfUPPRBVfSl , ST , $RH 

REAL K ,KAK ,L AM BO A ,LHA X , MH , MLE , MR , MSL , MSP , M V , MV I Ml 
A12= 2-/HI 1)/Hi 2) 

A34= 2 ./Hi 3 J /HI 4 ) 

AZ = A12+A34 

B12= i H12I-R* 1) >/RZ/IHIlH-H 121 t 

834= 4 81 4)*rU 4 *-*< 3l*tW3) ) /BZ/RZ/iH |3J+H(4)I -$ A L C I MJ/KM4 I M ) 

A( IP , 14 = ( Z./H4 U+B12 j/AZ/<H< U+H42) I 
AC IP, 2) - A12/AZ-A4 IP , 1) 

At IP, 31 = i 2 • /HI 3 ) +0 34} /AZ/(HI3)+HI4) I 
At IP, 4) = A34/AZ-AC IP ,3) 

M IP! = -TWhl*bZ*KZ^SALU MJ/AZ 

RETURN 

END 


SUBRUUHNE BURY12 t I , IM,I T,SURF t $URF&Vh 
C 

L B0KYI2 CORRECTS VALUED COMPUTED BY HK8 WHEN A VERTICAL MESH LINE 
C INTERSECTS A BLADE 
C 

COMMON /CALL UN/ MB I M 1 , MB I P 1 , MB QMl ,HBUPl ,M012M1 , MBl 2 PI , Md02 Ml , 
IMBU 2P 1,MMM 1, HM 1 , HM2 , hM 3, NT ,L)TLR, OMlK , PI TCH ,C P , EX PON ,T WW , C PT IP, 
21 GROG, TO I , ToC, LAMBDA ,TWL,I TMl N , I IMA X , NI P , I MS i 4 1 , BV |4| ,MVC 103 K 
31 V( 10 IK 1 TV! IOC, 4) , T Vi 100,41 ,0 TO MV iloO ,4) , BET A V i 13 0 , 4 ) , 

4MHC 10 0 t 4 ) , 0T UMH I iUQ,4) ,Bfc f A H i 1 00 ,4 1 ,RMHtlOO,4) ,8EHU03,4), 

SKM{ 10 0 KBfci 1U0) , DOOM t LOO I , SAL I 1001 , AAA UOOL 
CUMMGN /KHOS/ K HOHB l 1 GO, 4 I ,ftHO VB U DO ,4 1 

COMMON / HR BA AK / H C 4 ) , R 1 4 1 ,d 1 4 1 , KAK ( 4 1 , KA 44 I , I H i4 ) , RZ , BZ 
INTEGER B L U A T , A A N 0 K , t R SGR , S T RE N , SL CRD , SOK V L , A A T E MP ,SJRF,SUR F GV , 
IF IRST , UPPER, UPPftBV, Si , ST»SRW 

HEAL K ,KAK ,L AMGDA , LMA X ,MH , MLt , MR ,KSL , M SP , MV , MV 1 Ml 
Hi I J - ABSCFLUATU T ) T- T V( I M, SURF 1 > 

Rl I *= KHUVBi IM , SURF ) 

KAK I i ) -B VC SURFBVJ 

KA l 1 1=1 

RETURN 

END 
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SUBROUTINE »ua¥34U # IM*$URF * SUkFBV) 


C 

C B0RYi4 CORRECTS VALUES COMPUTED BY Hftb WHEN A HORIZONTAL MESH LINE 
C INTER $ELT S A BLADE 
C 

COMMON /CALCUN/ MaiHl, HdlPl ,MBOML ■M&QP1 , MB 12 Ml t MSI 1 Pi , M302MI , 
IMB02P IfMMMl* HMl # HM 2 * HH 3 * H T , 0 TL R ,DMLft t PITCH # C P , EX PQN iTtfW , CPT IP# 
21 GROG# TBI r Tb G r LAMBDA , ?RL # I TM1 N # I TMA X # N I P # i M$!4) , GV t4> # MV (100 It 
31 V f 10 1 )f I TVI IOC * 4) t TVi LOO f 41 ,L>TDMV(lQt) #4 } # BET AV (130 *4 J * 

4MHI 100*4) r U T DM K ( 1 0 Q , 4 ) ,6£TAHU00 #4 } *KMHU0Ot4) T BfcH 1133,4)* 

5RM4 100 USE* 100 ) t Qb U M ( i 00 i .SALtlOOl #AAAUt>0) 

COMMON /RKOS/ KHQHB U 00* 4 ) P RHUVB U OQ ,4 J 

COMMON /HRBAAK/ H ( 41 1 R 14 1 i b 1 4 ) f KAM4) # KA 44 > # I H 14 I # RL # iiZ 
IN I EGER SLOATtAANOKf CR SUM# STRFNtSLCRO i $Uft V L*AA TE HP *SURF #$U RF BY * 
IF IRST* UPPER# UPPRbV#$l »ST, %RH 
REAL R #RAR # L AM 8 DA * LMA X P MH t ML£ j MR »MSt P M$F # M V , MV I Hi 
iHt SURPJ^IHI SURF i + I 
IHS^IHISURF) 

Hi i } =A BS1 MVt IM HMH( 1HS,SURF) I 
Ri l )*fthOHBi IHS, SURF) 

Bt 1 I HS P SURF I 

*AKt I ) *6¥S SURFS V) 

*At I J = 1 

RETURN 

bNQ 


SUBRU U 1 IN E SUM 

SDK SOLVES f HE SET OF SIMULTANEOUS EUUATEGNS FOR THE STREAM FUNCT I UN 
USING THE Me: THU 0 UF SUCCESSIVE U V£ K -ft E L A X AT ION 

COMMON /AUKR HU / A t 2000 #4 ) t U( 2000 ) # K i 2000 ) # RhtO I 2003 } 

COMMON /iWPf GAMiAR , T I P # RHUI P , WTFL * R TF L$ P , UMEO A , QRF , BET A I # BET AQ # 
iMtH * MbU# M B 12 f HBU 2 rMM t Ndb i # NB L # NR SP , MR I 5G I # RMS P 150) * EES P I 53 ) , 
2BLDAT# AAn JKj ER SUR# SfRPN # 5LCRU P INTVL,SURVL 
COMMON /CALCON/ m IM L P MB l P 1 t m UMi # MBOP 1 # MB 1 2 Ml * MB! 2 PL , MBG2M1 , 
1MBG2P 1,MMH1# HM1 ( HM 2 ,HM 3 * M T *D TLR , DMLR ,PI TCH ,C P # E X PON ,T MU , CPT I P# 
2TOROG# IBI t TB 0# LAMBDA T T WL , I TMI N p l TMA X #NIP t I M5I4 ) *BVC4I P MVilD3 ) t 
iIV{ 10 i U 1 TVi 100 # 4 ) *TV1 100*4) #0 TO MV ( 1 00 ,4 ) f uET AV ( 13 0 r 4 } , 

4MH( 10 0 j 4) P OTOMH{ lQU t 4) ,BtTAHUU0t4l f KMHU30»4) #8EH(lD3#4i# 

5kM l 100 \*&hi lQOJtUBDMl 100) # SA L 1 1 00) *AAA(100) 

INTEGER BU>AT ( AANOR,fcRSUR f STRFN, SLCRD # SUR V L # AAT EHP # S J RF* SU RFBV * 

IF 1R ST f UPPER* UPPR8V # SL r ST, SRW 

REAL R ,RAR#LAHBUA,LMAX f MH , MLE , MR ,MSL ,H SP , M V , MV I Hi 

AATtMP = AAV UR 

IF t UR E .GE wd*\ uKF “0* 

IF IUKF,wT,l.) GU TO 20 
UR F ^ 1. 

UR FUP T = 2* 

10 URFTEM -UR FUP T 
LMAX = Oi 

20 IF i AATEMP *uT.U) HR l IE ( 6 t 1 Ul 0) 

ERROR = 0- 

SOLVE MATRIX EQUATION SY SOR # UR CALCULATE OPTIMUM OV ERRE LAX AT I ON 
FACTOR 

IP = 0 

DU luO i H ' 1 , MM 
IPU * XV £ IM) 

IPL - IV < 1M+1I-1 
IT = I TVi lM f 1) 

IF IIM.GT-MBUi IJ-iTVt 1M#3) 

IF { A A TEMP *G f .0 ) HR i TE 16 , 1020) IM#IT 
u-o 100 IP-IPU, EPL 

i Ft I T . GT . I T V { i M p 4 ) * A ND * l T* LE * I T V U M # 3 i 1 IT = lT+ITVilM,3l 
1- IT V( IM# 4 1~ i 
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IP l » IP-1 
IP2 * IP + 1 

C CORRECT 1P1 AND iPZ ALUNG PERIODIC BOUNDARIES 

I F ( IM.GE.MBi .AND * { M. LE .MB02. AND* U M. LE * MBO. OR * I M.G E .MB 1 2 M GO TO 33 
if( iT,60,irvii« t i)*0R.ir.eg*irvii«t3n epi = ipi+nbbi 

I FI IT . EQ . I TV { IM , 2 ! *Oft » 1 T* E U. I T V ( 1 M ,4 J I I P2 ^ |P2~NBfH 
30 IT3 ^ |T 
IT4 = IT 

C CORRECT IT3 AND IT4 ALONG LINES C-D AND K-L 
I Ft IN . NE.MBOP 1 J GO TO 40 
I Ft IT *LT * ITV t 1H- 1*111 IT3 = 1 T*N6UI 
40 1F( IH .NE.MBO J GQ TO 70 
IFtMB I 2-M8U J 50, 5G*60 
SO tF( IT .G£* ITVt IH,3)J IT4 ~ 1 T~NBBl 
GO TO 70 

60 If I IT.gT.ITVT IM*l,4l J IT4 - I T-NSBl 
7 0 IP 3 = IPF( IM-1,IT3J 
IP4 = IPFt IM + I, IT4) 

IF iOftF.GT.UJ GO TO 80 

C CALCULATE NEW ESTIMATE FOR LMAX 

UNEW = At IP, U*U( IP1J+AI IP ■ 2 J *Ut IP2 UA(IPi3)*UiIP3)+AI IP, 4 J*Ut I P4J 
IF ( UN i=W .L T ■ 1*£~2 5) UUPJ ^0. 

IF tUt IPUEU.O. J GO TO 90 
RATIO = UNEW /U l IP I 
LM AX- ANA X It RA T 1 0 , LHA X I 
Ut IP J = UNfcW 
GO TO m 

L CALCULATE NEW ESTIMATE FOR STREAM FUNCTION 8V SOR 

80 CHANGE = GRF*iKUP)-U(lPH-AUP,il*UllPlH-AUPt2>*Ut IP2 )+AUP*3|* 
iUt IP 3MAI IP* 4 J *U( I P 4 1 I 
ERROR = AM AX 1 1 ERROR ,ABSfCHANGE I ) 

Ut IP I - ut ip h-change 
90 IFl AAT EMP .Lfc *0J GO TO iUO 

WRITE (6, 1030 J II,lP,iPl ,1P2,IP3 rIP4 ( tAIJP,IUI = l,4l ,KUPJ 
10 O IT » I T* A 
AATEMP = 0 

IftOKF *GT ,i. > GU TO 110 
UR FOP T * 2,/t 1* +SQR T (AB SI 1 * — LMA XI J) 

H ft I T £ t t r 1040 I URFOPT 

IFt URF7tH~URFDPT*GT..0QO01.Uft, QRFUf T.GT, 1.999J GO TO 10 
WR LTE 16* 1O0OI 
ORF = CKFOPT 
GU TO 20 

UO IFIERSOR.GT.OJ Wft 1 TE ( 6 * 1 OS 01 ERROR 
I Ft ERROR, G I, *0000011 GO TO 20 
IN STRFN.Lfc.OJ KETURN 

PRINT STREAM FUNCTION VALUES FUR THIS ITERATION 

WRITE i 6, 1060 J 
MBIT = M IN Qt M &U , MB I 2M 1 J 
DU 120 IN =l,MbIT 
IPU = iVt IMJ 
IPL * IV( IM*1J-1 
IT VU - ITVUMtl) 

WRITE 1 6, 1020 lM t |TVU 
120 WRITE t 6, IU7 U J I U t I P I , I P = I PU , I PL J 
IFIMBI2.GT.HB0 J GU TO 140 
DO 130 lH=NtU2,MBU 
IPU ^ IVi EHl 

IPL - IVf IM) + ITVUH'4J-£TVC1M,11 
IT VU * ITVUM* 1) 

WRITE 16,1020) Irt.iTVU 
WHITE {6* 10701 tUUPUlP^I PUtIPLJ 
IPU * IPL + 1 
IPL ^ IV< IM*U-1 
f FI IPU.GT.IPLJ GO TO 130 
I T VU = ITVi IM* 3 ) 

WRITE «6, 1020* EH, I T VU 
WRITE t 6, 10? G J iUI IPJ, IP = IPU,IPLJ 
130 CONTINUE 
140 UO 150 IM =MBOP 1 ,HM 
IPU * IV( IMJ 
IPL = iVllM+l)~l 
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tfv u ^ irv(iM,3) 

WRITE i 6, 1020) IH,ITVU 
150 WRITE (6,I07G) ( U ( IP 1 , l P =1 PU a PU 
RETURN 

1000 FORMAT I 1H1) 

1010 FORMAT i a ani. IT IP IPl IP2 IP3 IP4 Atu A(21 

1 Ai 3 ) A(4) Ki 

1020 FURHATI5H IM *♦ 1 4, 6X f iiHl Ti * ,14* 

1030 FORMAT! IX, 14,516, 5F 10, 5) 

1040 FORMAT { 24H ESTIMATED OPTIMUM ORF = *F9,6I 
1050 FORMATION ERROR = ,FU,8) 

1060 FORMAT ( 1H 1 , 10A , 22H STREAM FUNCTION VALUES) 

1070 FORMAT i2X,lQF13*a» 

END 


SUBROUTINE SLAX 

Si AX CALLS SUBROUTINES TO CALCULATE RHO*W-SUB-M THROUGHOUT THE REGION 
AND ON THE BLADE SURFACES, AMO TO CALCULATE AND PLOT THE 
STREAML IN £ LOCATIONS 

COMMON /AUKRHU/ A ( 2000 ,4) * Ui 2000 ) , K ( 200Q I , RH 0 12003 ) 

COMMON / INP / GAM , Aft , T I P , ft HU I P , WfFL , n#TF LSP , OMEG A , (JRP , DEI AI f BET AO, 
lM0ltMau,Mbl2,MaO2,MMrNB8l , NB L *NR SP , MR * 50 1 ,RMSPC50) .BESPiSO )t 
2tiLDAT, AANDK, £R SUR , Sift FN , SLCRD * i N TVL , 5URVL 
COMMON /CALCO H/ Mb 1H 1 , MB I P L , MB 0M1 , MB0P1 , MB12 HI , MBI 2 PI , MBU2M1 , 

1HBU 2P itMMMl, hMl,HM2,NM3f HT »U TLft , DMLft ,PI TCH ,C P , EX PUN ,T W W , CPT IP, 
2TGRUO, Tbl , Tb O , LAMBDA , TWLtI TMI N , I TMA X ,NI P , I MS (4 1 , BV i 4 1 , MV ( 100 I , 

31V i 101 I, I TV! 100, 41 , 1 VI 100,41 ,0 TO MV (100 ,41 , BET AV 1133 ,4) , 

AM Hi 100,41, DTUMHt iUO, 41 , BE TAN U 00 ,41 ,RMH (100,41 ,B£H(UD ,4) , 

SftM( 100 J,bEU00l ,PBOM ( 1 00 1 ,SAUiOQ) ,AAA (1001 
COMMON /SLA/ T SU BOO 1 , Ui NT ( 81 
DIMENS ION HSLi 100) ,KKM18I ,P(41 

DIMENSION w( 20001 ,RWM( 20001 ,BETA(2O0O) ,WM8U3Q ,41 , WT Bi ID 0 , 4) , 
1XD0WNC £00 1 , Y AC RO S ( 800 1 

LOUlVALENCt (Ai 1 , 1 J , W ( ill , 1 A ( 1 ,21 , ft W M ( l ) 1 , ( A ( 1 ,3) , BET Ail) )» 

XI A( 1, 4),hHB| 1 ) 1 , ( A ( 401 ,4) * WTB ( 1 ) ) , (A (SOI *41 ,XDGWN(lf 1, 

2(M l),¥ACKOSCiH 

INTEGER B L DA T , A A MDK , E ft SGR, STRFN , SLCfU) , SURVL, AATEMP ,S JRF ,S URF BV , 

IF Ift ST, UPPER, UPPRb V , 51 1 ST, 5RW 

REAL K ,KAK ,L AMBOA , LMA X ,MH , MLE , MR , MSL *H5P , MV, MV I Mi 
DATA ( KKX ( J ) ,0^4, 1 B * 2 ) / 8* l H * / 

CALL 5LAVP AND SLAVaS THROUGHPUT THE REGION 

irvu^ irvtiti) 

ITVL^ IT V ( i, 21 
DO 10 iMslsMalMi 
10 CALL SLAVP (IM, I TVU, ITVL J 
HBOT* MINOiMbO ,HBI 2MI) 

IF (Mb 1.GUMBUT) GO TO 30 
DO 20 IM-Mb I ,MbO T 
1= 0 

20 CALL SLAV BBt IM, 1,2, 1,1 ) 

30 IF (MBUP l.GT .MB I2M11 GO TO 50 
ITVU* IT V ( MBUP 1,3) 

ITVL= IT V ( MBOP 1 ,4) 

DO 40 IM 80 P l , M B I 2M 1 
40 CALL SLAVPUM, ITVU, I TVL) 

50 IF (MS I2.GT,MbUl GU TO 7D 
DO 60 IM-Mb I 2,MB0 
1= 0 

CALL SLAVBQt IM, 1,4, 1,1 I 
60 CALL SLAVBUI IN, 3, 2, 4,1 ) 

70 MBGT= MA X 0 ( M BOP 1 , MB 121 

IF 4 MBOT * GT *MSQ 2 ) GO TU 90 
DP 80 IM= MB OF ,MBQ2 
1^ O 
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60 CAUL SLAVBBI IM , 3, 4* 3 , i I 
90 ITVU* ITv(MBU2PI»3l 
ITVi* 1 rviMaU2Pl,4| 

DO 100 IN^H602P1,MM 
100 CALL SLAVPUM.lTVUtlTVLi 

PLOT STREAMLINES 

IF CSLCRD.LE.Ol HE TURN 
DO U0 IM^UMK 
110 HSLI mi - HV( IM ) 

KRKt 1 1 « 7 
KK*U 2 > - B 
HKKI3I * HH 

pi n - i. 

PI 31 s 0* 

PI 4 I = 0 * 

WRITE! E, 1000 1 

CALL PLUTHYI MSL »TSLtKKK,P| 

WRITE! £, 1010 1 
RETURN 

1000 FORMAT 4 2HP F 1 3QX ■ 1 &HSTREAHH ME PLUTS 

1010 FORMAT 12HPL , 40X, 7GH STKE AML I N£ S ARE PLOTTED WITH THETA ACROSS THE 
IP AGE A m M DO*N THE PAGE) 

END 


SUBROUTINE SLAV 

SLAV CALCULATES KBU* W-SUB-M THROUGHOUT THE REGION AND ON THE BLADE 

SURFACES g AND CALCULATES THt STREAMLINE LOCATIONS 

COMMON SR W, I T tR 1 1 END • L EK ( 2 1 ,NERf21 

Cummin /AUHRHj/ A 4 2000 ,41 , Ul 20001 ,K(200O) ,RHO(2D03 1 
COMMON /IMP/ GAM,AR , T IP»RHOI P ,WTFL , WTP LSP , QMEG A , OR F, BET A 1 , SET AG, 
iMB I,M6Q,M8I2,N0G2,HH P NBBI , MB L ■ NR $P , MR ! 50 1 , RMS P 150 1 ,6ES PI 53 l , 

2BL OAT * AAN DK t EkSuK,$TRFN, SlCRD * I N TV L * SURVL 
CUMMIN /GALLON/ MB i M 1 , MB I P 1 , MS GM1 , MBQP1 , MB 1 2 Ml , MB 1 2 Pi , MB02 MI * 

IM BU 2P i ,MMM 1, HfU, HK2 ) HM3»htT,0TL« ,DMLK ,PITCH ,C P t EXFON,TMW , CPT I F, 
2TGRUG, Ifli ,T6u, LAMBDA , ! WL , I TM I N , I THA X, NIP ,1 MS 14 1 ,BV44I ,MVUO0 )t 
AlVi 101 >,ITV( iOC,4t,TV( 100,41 ,3 ID MV 4 1 U0 ,4 1 , BET AV 1 130 ,4 I , 

4MH( iU0,4l ,DTDMH( 100,41 »6E TAH 4 100 ,4 I ,RMH41D0,4l ,BEH!10 3 ,4 I, 

SRM 4 IDO I , B 1 1 1 00 ) , Od 0 M ( 1001 , SA L 4 1 Q01 , A AA (1001 
COMMON /SLA/ T SL 4 BOO 1 * ill N T 4 8 I 
D I MENS ION TSPl SOI ,USP 4 SOI ,0001(301 , Ti NT 481 

OIMENS Ion Hi 2000 J,RWM I 20001 ,8£ TA C2G001 ,WMB U0O ,41 , WT B 4 13 3 , 4 1 , 

LX DOWN ( 600 l.YACROSi 8001 

EQUIVALENCE 4A4lfll,W(lJ 1 , I A U ,2 1 * RUM 1 11 1 , t All ,3 1, BET Mill, 

UAL 1,41, WM64 11 I r I A I 401,41 ,KTBUJ ) , (A 4801, 41 ,XDUWNUI I , 

2! Hi it, VACH0S4 1) 1 

SLAVP CALCULATES ALONG VERTICAL MESH LINES WHICH 00 NOT 

Intersect blades 

ENTRY SLA VPi 1M , I T VU, I T VL 1 
LUC 35 0 
il = 0 
Ni- 8 

N5P- I TVL - IT VU+2 
IP - tvdrt 1-1 
DU 10 IT* l, MSP 
IP * IP + 1 

fSPUTl * FLOAT! 11 + 1 TVU-1 1*HT 

10 USPI IT 1- UIIP1 

USPfNSPl * ySP m + i* 

IP = I Vi i M I 

INTO * IN T I U C I P 1*3*1 

IF IU4iPI.GT.un iNTU=lNTU+i 
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00 J = 1p 5 

UiNltJ) = FLOAT! INTU) /3, 

20 INTO - INTU+1 
UINT 16 i = B VI A* 

30 If IU INTI 6KGE.UUP H GU 10 40 
U IN* T I 6 )= UINTI 6 i * i * 

GO TU *0 

40 If OilNTI 6l.Lfc.UliP Hl.J GO TO 50 
UlNTi6>= U1NTI6>-1. 

GJ TO AO 

50 OINK 7 UIN H U 
OINK S )= UInTI 11 
GU TU JOO 

SLAV Bd CALCULATES ALONG VERTICAL ME SH LI Nfc S WHICH INTERSECT BLADES 
ENTRY SLAVBBi Irt, UP PER t LOWER t UPPftBV , 1 1 

INTEGER t>LOAT T AAND&,fcftSGR* STRf N # SLCRO ,SURVL, AATEMP rSURF ,$DRPSV i 
IF IKS T, UPf EH, UPPRBV, SI * ST t SftW 
REAL K t KAR * L AM BOA » LHA X rMH> i MLE * MR , N SL t H 5P , H V * HV I Ml 
LUC* 1 

1TVUP1 * I TV i iM f UPPER i 
ITVLM1 * I TV I IM # LOWER I 
IT VU * ITVUP1-1 
IT VL ^ iTVLMl+1 
NSP = ITVLMTVU+1 
TSPt 1 I = TV! IM, UPPER 1 
TSPtNSP 1 = TVt 1M , L UWER 1 
USP( U = BVI UPPRBV) 

USP(NSP) = d Vi LOWER I 
IP = IPF« IM t ITVUPU-1 
N$PN1 * NSP- I 

1 Ft 2,Gf .NSPM n GU TO 70 
DU 60 IT^2*N$PMI 

IP = IP*1 

TSPtin = FLOATU T+I TVU-l |*HT 
60 USPi IT* = 04 IP 1 
70 11 = I 
1= U + J 

UINTI 1 )= 6 Wi UPPftBV* 

INTO * I N T £ U IN T I I 1*5, I 
If (U1NT4 D.Gt.O.i I N T U-I N TU*l 
SC l * I+I 

UINTI l i * FLJA ft INTU1 /5. 

INTO = IN T U* l 

IF 10 INK i KLT,BV<LUWEK> J GO TO 00 
U IN ft I 1= BV£ LOWER ) 
if iLOWEk-UPPEK *NE.li GO TO 90 
I = 7 

UINTI l i = BV 4 A i 
90 U IN It S I * SV t LOWER 1 
N I* l- II 

IF INI *EQ * 7} N I * 0 

FOR BOTH SLA VP AND SLAVSB , CALCULATE RHQ*W-SUB-M IN THE REGION, AND 
RHJ*W AT VERTICAL MESH LINE i NFE RSE CTI 0N$ ON THE BLADE SURFACES 

ICO CALL 5PL INC ( TSP , USP *N SP #0 UO T * AAA J 
FIRST* t l-LUC!*ITVU+UJC*l TVUPI 
LAST - 4 l-LOCi*lTVL+LUC* I TVLM1 
IFI FIR ST. GT, LAST) GO TO 120 
IT * LUC 

IPU * iPFt IM* FIRST > 

IPL * IPFt IN, LASH 
00 110 IP * IP U? IPL 
IT = IT*! 

no RwMtiPt = ouyiti ii*wtfl/be t iM» /mu m 

120 IF tLQC.tU.O) GO TO 130 

WMBI 1M,UPPER I = OUDTi 1 i * WTf L /BE U Mi /RM tt Ml 
WMBt IM rLDWER J = DUDT4NSP **WTPL/BE£ I M * / RM t I Mf 
RMDTU2 = t RM t 1M )*Q TO MV ( 1 M , UP PE ft) ***2 
RMUTt 2 = (KHt lMf*DTDMVl IM, LOWER) 1**2 
If ( RM DTU Z.G T * 10000, 1 WMB(IM,UPPEftl * 0* 

IF IRMOTl 2.GT. 10000*1 WHa 1 1 M , LOWER) * 0. 
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non non nnnno 


MMtH IM, UPPER 1 = AbSflrfrttmH,UFPEftH *SaKTU .+RMDTU21 
tfMd( IM ,LJMEft 1 = AB S tWHdl I M , LUWER) * *SORTU * *RMDTL2 1 
130 If ISLCRD.Lc.01 RETURN 

11 = 11 * 1 

CALL SPL INK US*' ,TSP,NSP,UINT (III ,NI .TI NTUll f AAAI 

12 = N i*i 1-1 
Oil WO J*II, ll 
L- U- 1HMM* IM 

140 TSLtt 1 = UNTfJi 

IF (UPPER ,EJ .WANO.LJtfER.EU* 41 K£ I URN 
IF (1M.EU.LI rift I TE ( 6 , iOU D ) 

WR 1 T t < 6, ID 10 1 MV I IM1, ( UINTt Jl , TINT Ul , J=1 ,dl 
RETURN 

1000 FORMAT* lHl////30Xf22HSTKEAKU ME GOOROI NAT £ 5/ / / / 5 X, 12HM COORDINATE, 
13* 7 x, i onstream fn* ,i ox, shine ta ,4xi / /i 
1010 FORMAT* IX, TGid.TZl L9X, oGId.TH 
END 


SUBROUTINE r ANG 

TANG CALCULATES ft HU* H- SUB - THE FA A NO THEM RHU#tf THROUGHOUT THE REGION 

AND UN THE Blade SURFACES, AND CALCULATES THE VELOCITY ANGLE, BETA, 

THROUGHOUT THE REGION 

COMMON SRwt I TER , I END , L LR I Z 1 ,NErt (21 

Common /auRRhj/ a * 20uo ,41 *u(200u wk( 20G0) ,R4a(zJ03 \ 

COMMON /INF/ GAM , A ft * T I P , RHl) I P , M TFL T InlTF L5P , UME G A , GRF , BET A 1 , BET AO , 
IMdl, HBO* MB 12 ,HQ(J 2 * MM , NBb I , NdL , NR SP * MR ( 50 I , RMS P 150 J , BESPf 50 I, 
20LOAT, AANOK, £ft SGR, STKF N , SLCRU , I N T V L , SURVL 
COMMON /L ALC UN f Hd IM1, MB I P l , MB0M1 . MdOPl ,M8I2M1 , MBI 2 P 1, M8Q2 Ml , 
iMb0 2P 2, HM 3, NT, DTL ft , PI TCH ,£ P f EX PUN t T tfW , CPT 1 P, 

2T G ft O G t IB I , TB (J , L AM BOA , T WL , I 7 MI N , I IMA X ,NiP ,1 MS !4 1 , BV (41 , MVUOO }, 

3 IV { 101 I, ITVt 10 0 » 4 W TV ( lOO ,41 ,D TD MV ( LOO ,4 1 , BET A V ( 1 T O , 4 I , 

4MM 100,41 ,OTUMH( 1U0,4 I ,tttTAHU0O ,4 1 p KNHfiao t 4l ,0EH( 103 ,4 ), 

SRH { 100 I,d El iOOJ ,DBDM( 1001 ,SALllOQ) ,AAA U0G1 
D IM EN S ION SP M 1 1001 ,USP (1001 ,DUDM(lUGl 

DIMENS ION til 20UQ* ,RtfM( 20Q0I ,BE TA (2 0001 i wMB i 1 00 ,41 ,tfT6I 133, 41 , 
IXDUtffcM BOO 1 , VACRUSt dOOl 

Egg (VALENCE 1 A t 1 , 1 1 ,tt ( U I , 1 A i 1 ,2 1 , Rtf Ml 1 1 , WU ,3 I , BET A ( I 1 1 , 

it At 1, 4 1, tfMB ( 11 I , (A [ 401 ,41 tHTBUl 1 , 1A ISOl ,41 ,XDOHN( 111, 

2(R( 11, YACR0S4 1) 1 

INTEGtK d L DA T , A A N U K , b K SDK , STRF N , SLCRD , SUR V L , A A TEMP ,S J RF ,SORFBV , 
IFIRST , UPPER, UPPRBVtSi ,ST,SRrt 
REAL ft ,KA r,L AMBDA,LMAX,MH,MLE , MR ,MSL , MSP , MV , MV I Mi 
LOGICAL A DDL , ADO 
EXTERNAL BL1 ,dL2,BL3,8L4 

PERFORM CALLULAT IONS A LUNG ONE HORIZONTAL LINE AT A TIME 
IT = IlMiN 

1U IF (IT *GT « I TMA X 1 RETURN 
AUDL = .FALSE. 

ADO - .FALSE* 

III = IT 
SI = 0 

ON THE GIVEN HsJKlZONTAL MESH LINE, FIND A FIRST POINT IN THE REGION 

IF! I T * GE * U *A ND • I T • L T . N b B I ) GO TO 60 
IM = MBIM 1 

20 IM = IMfl 

IT V 1 =* IT V4MdO,3J-NBdI 
iF(Mdl2*GT,Mb01 I T VI - 1 TV ( MB OPl ,3 1 

IF( IM .EO.MbJ .ANO*I T*G€*I TV1. AND. IT, LE. I TV (MB0,2I -NBtU 1 GO TO 2D 0 
IF! IM .GT.MBO 2P11 GU TO 200 
DO 40 SUR f = 1 , A , 2 

If i !M *GT *MdO. AND. SURF *EQ. 11 GU TO 40 
I TV IM I * 1 TV ( iH- 1 , SURF 1 
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1 H SUK f *£ D .3 .AND* I M.EQ .MB0P1 } ITVIM1 * ITVIM1-NBBI 
IH IT i *UE .IT V(IH, SURF! .AND. I TI.LT.l TVI Mil GO TO 7 0 
40 CONTINUE 
SURF = 1 

IH IH * FLf .M bUP 1 .AND . I T * EQ . I T V ( Mb U ,1 ) -1 , AND. I T V ( MBU, 1 1-ITV(M6G,2 } 
l+NBBi.EQ.21 GD TU 70 
00 50 SUR F = 2 t 4# 2 

IF i IM.Lt -MB 1 2 .AND * SURF - £ Q . 4) GU 10 50 

IF 4 II .L t * IT V4 IM , SURF J .ANQ.1T.GT.I TV( l M-I , SURF J 1 GO TO 70 
50 CONT IN L£ 

GO TO 20 
C 

C FIRST POINT IS ON BOUNDARY A- N 
C 

6Q IH I- l 
IM = 1 

SPHI 1 1 = MVt 11 
USP (11= Ui IT+iJ 
GU TO ED 

c 

C FIRST POINT IS ON A BLADE SURFACE 
C 

7 0S! = SURF 
IT 1 = IT 
ADO - .FALSE . 

IF I AODL . AND. SI. EU. 31 ADD - .TRUE* 

IF i AUDI ITI^I T-N8BI 
IM1 = IM-1 
IN 2 - IH 

TH = FLOAT! I TI 1*HF 

IF ( SI .Eu.i.AND* IMl.LT.M&UJ TH = TH— FLOAT 4 NBB 1 1 *HT 
MV IH1 - HVt IMH 

IF 4 1 IM.tU.MBIP l.AND. { SURF . EQ, 1 . OR . SURF. E 0.2 1 1 . OR. ( I M. 6 Q * MB I 2P 1 
1 .AND* I SUR F. £0. 3, OR, SURF. EU. 4 11 1 MVI Ml = MVI Ml +■ I MV (I M2 1- MV I Ml 171300* 
L ER4 21 “LO 

€ BLCD (VIA RUuTI CALL NO. lo 

IF ( SI *EU . 1 . AND * IM 1 *NE . MB U J CALL ROOT! MVI Ml * MV ( IM2 1 ,TH , BL1 , OIL R, 
iANStAA A1 
LERI 2 1=11 

C BL CD (VIA ROOT 1 CALL NO# U 

IF ( S 1 *tu . 3. AND • IM l. NE -MB02 1 CALL ROOT I MVI Ml r MV ( I M2 1.T H T BL3 , DT L R , 
IANS* AAAI 
LERI 2 1 =12 

C BL CD (VIA R3UI1 CALL NO. 12 

IF ( S 1 .EU* » 2 ) CALL ROOT (MVI Ml, MVI I M2 1 ,TH , B L2 , DT LR , A NS * AA A 1 
LERt 2) =13 

C BL CD (VIA RJQTJ CALL NO. 13 

IF4SI.E0.A1 CALL RUOFt M VI M l , M V I I M2 J t TH ,BL4 ,QT LR , AN5 , A A A ) 

IF( SI • EQ . 1 .AND. JM1.EQ.H3U1 ANS = MVIMBGI 
IF( SI .EQ.3.ANU .IM1.EU.MBG21 ANS - MVIMB02J 
SPMI IM 11 = ANS 
USP 4 iM 11= 8V I SI 1 
C 

C MOVE along horizontal mesh line until mesh line intersects boundary 

c 

B Q IT 1= IT 

90 I TV 2 = IT VIMbO ,3) 

IF (MB IZ.GT.MBO} I TV2 = I T V I MB OP 1 ,3 1* NBBI 

IF 4 IM.Ew .MB UP l.AND. I T.GE. 1 TV2.ANQ. I T. LE.I TV (MBU ,21 1 ADD- .TRUE* 

IF 4 AD L 1 A DDL = .TKUt* 

IF ( ADC1 I T 1 =1 T-NBBI 

IF 4 AUL.ANU. Sl.EQ. 31 USP U Ml J =B V * Si) *1 * 

IF 4 IN .LT .MB I.OR .IM.GT.M3U21 SO TU 120 
DO 100 SUR F = 1*3.2 

IF I IM *L t .MB 12 .AND. SURF. £ w *3 J GO TU 100 

IF 4 I M *E j . IM 2 .AND. Si * E 0 . 4 . AND. SURF * £ Q. 3 1 GU TO 103 

I T V IM 1 = I TV I IM-1, SU*F 1 

If i IM .eg .MBUP l.AND. ADO 1 ITVIM1 » ITVIMI-NBBI 
IF 4 IT I*L T .1 TVI IM.SURF 1.ANU.I TI.GE.i TVIHU GO TO 140 
100 CONTINUE 
SURF * 3 

IH IM .EU.MBI 2.AND, IT.EO.I TVI MB 12 ,31-1. AND. I TV < MB 12 , 3 I -I T V ( MB 12 , 4 I 
1.EQ.21 Gi TO 140 
DO 110 SURF = 2* 4 , 2 


ii= f fm*gt .MB 0 *ANQ*simF*ec* 2 ) go ro no 12 s 

If Uf^tG* TM2* AND*Si *tG* 3*AMD-SURF -EQ-4J GO TO 110 ! 2b 

1 I V 1 Ml = ITVUM-l,Sl>RFl 127 

IF UM-EQ-MfcUPUAhD. ADD* ITVFM1 « 1TVIM1-NB&I 128 

IF UT I *GT*1 rvtlM, SURF ) -AND, IT 1 .LE . I TVIM I ) GO TO 140 129 

110 CONTINUE 130 

120 SPMI IH) = PVi FM1 L3! 

IP - IPR IK, ITU 132 

USPtmi = UtlPl 133 

IF IAOO) USPIIM1 - USPUM1+1* 134 

IF UM*EQ.MM) GO TO 130 133 

lit* w + l 136 

GO TO 90 137 

C 1 36 

C FINAL POINT IS 0* BOUNDARY C-H 139 

C 140 

130 IMT - MM 141 

GO 10 bO 142 

C 143 

C FINAL POINT IS ON A GLADE SURFACE 144 

C 146 

140 ST = SURF 146 

I MTs fM 147 

|Mm= lMf-1 148 

TH - FLOAT UTI I*HT 149 

IF F$T,tQ,3,AND,iHr.Lt.MD0) TH - TH— FLOAT! NB Si )*HT 150 

MVlM * NVUMMU 151 

IF f < JMTMUEC.PUIJ.ANC,! ST*CQ.l-OR*ST-EO-2n 152 

l KVIM1 - MV[Ml+!M\H IMn-MVIMimOOD* 153 

IF F UMFMl*tQ*MBJ 2I,AN0* F$F, F0,3,QR*ST*Ea,4T -AND. 154 

1 I 1TV1MU12.31-ITVC M8 12,4 I * FO. 2* OR, I T V ( MB I 2 t4 J - I T Vf MB l 2 , 3) - EG* 155 

2 ^BU-2M HvIMl - MVlMl+tMVt fMT)~MVlMl)/lOOQ* 156 

LfcR F 21 = 1 4 157 

C BLCI) IVIA -TOOT! CALL NO* 14 l 58 

IF { sT. FN* 1 . ANU, I MT+VF * M8 J I CALL ROOT ( MV 1 M I , MV E 1MT J , TH*BL1 , 159 

IDTLR.ANSt AAA) 160 

LER 121 =16 161 

C MLCO {VIA RUUTl CALL NO, 15 162 

IF t ST.FO.S.ANu* [WT*NF *MB F2 1C ALL ROUT i MV I M l , MV l 1MT ) , TH, KL3* 163 

IDTLRtANSt AA a) 164 

LERI2J-16 166 

C GLCD (VIA ROOT 1 CALL NO* 16 166 

IF { ST-E0-2JCALL ROOT ( MVt HI t MVUMT ) , TH, BL 2, DTLR t AMS » A AA 1 167 

LEK12J = L ? 166 

C BLCD t V 1 A RUUT 1 CALL NO* IT 169 

IF t ST* £0*4 1C ALL ROOT I MV I Ml ,MVF I Mf 1 # TH*6L4*DTLK t AN5 , AA A } 1 70 

1 F 1 sT-FQ, L.ANU*inT*EO*MBn AMS = MV { MB I I 171 

fFt ST-FD, 3.AN0* IFT-E0-H612 1 ANS = MV ! MB I 2 ) 172 

SPMIMT) = ANS 173 

USPt !M11= □ V ( 5 T ) 134 

U !5T*EU-3.ANO.{MT*LE*MBU) USPt IMT > « 6V|4l 1T5 

[ f t UiD J OSPlIrtTl = USPtlFTJ + 1* 176 

IF !S1*FG*3* AND* ST ,E 0 . 2 ) USPUM11 = BVtSU + 1- 177 

IF ( S 1 ■ F 0 * 3 * AND *ST - E0-3I USP1IH1) = BVlSU+1, 170 

C 179 

C CALCULATE RHO*h-SUb-THFT A A NO THEN RHU*K AND BETA IN THE REGION 100 

C 1 81 

160 NSP* IMT-lMt+I 182 

CALL SPLINL { 6PM [ LM1) , USP ! IM1 i*NSP t OUDM{ iMfc 1 * AAA ! I Ml JI 1S3 

FIRST-1 184 

IF I {MI *NU * 1 1 F mST= 1 F2 1S5 

LAST- MM 166 

IF MMT-NE*KMi L A$T= 1 MTM l 167 

IF FRR$T,GI -LAST) GO TO 170 16® 

IT f = If 16*7 

IF {FIRST*GT-MB0P1.ANC*ALDI IT FM T-NBBl 150 

00 160 1=FIHST*LAST 191 

IF { ADD - AND - I - EG-MBQP 1 I ITI-IT-NflBl 192 

RWT - -DUOMUlHjTFL/BEUf 193 

IP = TPFU , ITU 194 

MIR - SORT !RWT**2+RV*M IP 1**2) 195 

160 mi* E IP) " ATAN2 tRHT, RHMIIP1 M67*295T?9 196 

C 197 

C CALCULATE RHO*N ON THE BL AO t SURFACES 190 

C 199 

170 IF | FM L -F0* 1 ) GU TO 1 BO 200 

CALL SEARCH ! SPM t I Ml) .Si , IMS I 201 

AMS ” -OUDMUwt L*WTFL/BEHUHS,SU 202 

ViThiiHSiSii - Aost ANS)*sofm i* + i*/tRMHi ms,si) *DTDKmiHS,sin**2) 203 

100 IFF l Ml *E0 *MM ] GO TO 2C0 204 

CALL SEARCH t S PM 1 MT J , ST * I HS > 205 

ANS = ”OUDM£ iHT)*wTFL/BEh{ FHS* ST J 206 

WT0UHS,ST} . auS C ANS F *S OR T t l * + 1 • / ( RMH ( I HS * ST 1 * D T|>MH ( 1 H$ t S T M * *2 ) 20? 

190 GO TO 20 208 

200 IT - IT+1 2®9 

GO TO 10 210 

END 211 


180 


214 

218 

222 

226 


25S 


279 

283 

286 

295 

301 

307 

313 


81 


SUBROUTINE SEARCH (DI ST, SURF ,1 S* 

C 

C SEARCH LOCATES THE POSITION OF A GIVEN VALUE OF M IN THE MH ARRAY 1 
C 

CO HHDH /C ALCUrj / MB IHi , M3 IP 1 , MBGM1 , HB DPI , MB 12 Ml T MSI 2 PI f MBG2 Ml , 

1HBO 2P l* MMM 1, Hrtl,HM2,HM3,HT,L>TLR,0HLR, PITCH ,£ P§ EX PON tTWW * CPT I P, 
2TGRUG . TBI,TBQ,LAMB0A,TWL,f TMIN,ITMAX,NiP.IHS«4l,BV(4l,MVliO0), 

3iVI 101), live 100,41, TV< 100,4) ,0TOMVUO0,4) , BET AV ( 13 □ ,4 I , 

4MHl 100.4*, 0T0MHU00, 41 ,BETAH(100 ,4J ,RMHUQ0,4* ,B£HU00,4», 

5RH( 100 l,BEl 1001 ,DBDMUl>0J,SAU10G> .AAA 41001 
INTEGER BL0A T . AANDK ,ER SOR , STRF N, SLCRO , SUR VL, AATEMP , SURF,SURF BV , 
if IRST t UPP ER, UPPR& V » SL , $t, srw 

REAL K ,KAK, LAMBDA, LrtAX,HH .HLE , MR ,M5L , MSP ,MV , HV I Ml 

00 10 1 * 1,100 

IF 1 AS SI M H( I , SURF l-DI S 11 , GT, DHLR 1 GO TO ID 
IS ^ I 
RETURN 
IQ CONTINUE 

WRITE (6, 10001 01 ST, SURF 
STOP 

1000 FORMAT I38HL SEARCH CANNOT FIND M I N THE MH ARRAY/ 7H DIST *,G14.6, 
I10X, 6HSUR F * , G 14, 61 
END 


SUBROUTINE VELOCY 
C 

C vaOCY CALLS SUBROUTINES TO CALCULATE DENSITIES AND VELOCITIES 
C THROUGHOUT THE REGION AND ON THE BLADE SURFACES, AND IT PLOTS 
C TFE SURFACE VELOCITIES 
C 

COMMON / A UKR HD / A 1 2000,41 . U( 2Q0QI , K (2 DOO ) , RH 0(2003 1 
COMMON / IHP / GAM , AR , T1 P ■ RHU1 P , WTFL , WTF LSP , OMEGA, DR F# BET AI , SET AO, 
LHG[,M8Q,H8I2,MB02,MM f N63I , NBL , NR SP t MR ( 50) ,RM5PB01 ,SESP150 l, 

2 BLOAT, AANUK.ERSOR, STRFN, SLCRD , I NTVL , SURVL 
COMMON /CALCON/ MB I Mi , MB 1 P I f MBOMl , MB OP l , MB 1 2 Ml , MB I 2 PI . MG02H1 , 
1H&U £P i ,MHH 1, HM 1 ,HM2,HM3,H T ,0 TLft ,DMLft ,PI TCH ,C P, EXPON ,TWM * CPT I P, 

2T GROG, mi, FBOtLAMBDA , TWL,I THI N, I IMA X, NIP ,t MS i 4 ) , S¥ 14 ) , MV ( 100 ) , 
3IV1 101 1, IT VI IDO, 4) ,TV| 100,4) ,0 TO MV (100 ,4 ) , BET AV I ID 0 ,4 I , 

4MH( 100,41 ,DTUMH( 100,4) ,8fi TAH ( 1 00 ,4 1 ,RHH (100,4) ,B£H (100,4), 

5RM( 10 0 ) , B £ f 10O),DB0M( 100) ,$ALUO0) ,AAAtIOO) 

DIMENSION KKKI18I 

DIMENS CON Ml 2000 1 . R WM ( 2000) , BE TA 120001 ,WMB (100,4) ,WT&(100,4) , 

IX DOWN I BOOT, YACR0S1 BOO) 

EQU I VALENCE (A i 1 , 1 1 , HI 1) I , ( A ( 1 ,2 ) *AWM(1) I , I A (1 ,3 * , BET A! 1) ) , 

II A* 1, 4i,WMBl 1I1,U! 401,4) tWTftfO ) « (A (SOI ,4) ,X0QWN(1 J i , 

21 K i 1), YACROSl 1) ) 

INTEGER BLOAT,AANOK,£R£Uft,STftFN, SLCRD , SURVL , A ATEMP , SURF ,SURF8V , 
IF IRST, UPPER, UPPR B V, SI, ST, SR W 

REAL R , R AX ,1 AM BOA . LMA X , MH , RLE , MR , MS L , MSP , MV, MV I Mi 
DATA KRK1 4J/1H*/,KKKI 61/LH0/,KKK(8)/1H=/ ,KKKI10)/1H(/, 

IKK M 12 )/lH+/,KKK( 14J/AHX/ f KKMl6) /lH*/,KKMiaj/lH)/ 

C 

C CALL VELP, VELBB, AND VELSUR THROUGHOUT THE REGION 
C 

CALL vaPU,M 8 IMi,l, 2 ) 

IF (NS I2.GT+MBU) GU TO 10 
CALL VtLBSfMeipHB 12M1, 1,2) 

CALL V EL &B 1 M SI 2 ,MB 0,1,4) 

CALL VEL8B(MB12,M80,3,21 
CALL V EL8&I M SOP 1 ,MBU2 ,3,4) 

GU XQ 20 

10 CALL V tL B fl 1 M 3 i , MBQ ,1,2) 

CALL V £L P 1 MBQP I , MB I 2M 1 ,3,4) 

CALL VELS0(N8I2iMBO2,3,4* 

20 CALL vaPI«B02P 1,MM,3,4| 

CALL VELSUR 

c 

C PREPARE INPUT ARRAYS FUR PLOT OF VELOCITIES 
C 

IF l SURVL *LE *0) RE TURN 
NP2 * 0 

C SURFACES 1 TO 4 - TANGENTIAL COMPOf€NTS 
DO 50 SURF “1,4 
NP1 = NF2 



n n n o 


IrtSS = IMS! SURF I 

IF I IMS 5 *L T * 1 i GU TO 40 

UO SO = IH SS 

IF I A&SiDTDMHi INS, SURF )*RMHUHSfSURF 1 1 .LT, .57735 1 GO T 0 30 
MR I * NPl + i 

YAtROStNPil * liTSI TH$ t SURF > 

XDUHN1NPU = MHi IHS, SURF 1 
31? CONTINUE 

40 KXK f 2* SUK F + l 1 - NPI-NP2 

50 NP2 = NR I 

C SURFACES 1 AND 2 - MERIDIONAL COMPONENTS 

m ao sur f *i , 2 

NP i = NP 2 

IF I MB IP l.GT .MBUM1 \ GO TO 70 
00 60 IH^MBIP LpMBUMi 

IF IABS( 1 >TQHV( lM,SJKJ-i*RMUttn. ST. 1.73211 GO TO 63 
NP1 = NP1 + 1 

YACUOSINPII = WMtHIM,SURFl 
XDUrtNl NPl } = MVt IMI 
60 CUNT tNUE 

70 XKK 4 2* SUR F *9 I * NP1-NP2 
BO NP2 * NPl 

C SURFACES 3 AND 4 - MERIDIONAL COMPONENTS 
DU 1 LU ' $URF = 3* 4 
HP l = NP 2 

IFIMBUP1.&T.M602MU GO TO 100 
DO 90 JM=MG1 2P1,MBU2M1 

IF t AB5( UTUHVUMtSURF i*RMUMl i.GT. i. 7321 1 GO T O 90 
NPl ^ HP 1*1 

YACKQSlNPLl = WMBUMpSURFI 
XDUWNt NP 1 1 = MVt IM J 


90 CONTINUE 

100 MtKi 2^ SURF + 9 } = NP 1-NP 2 
110 NP 2 = NP 1 


C 

C PLOT VELOCITIES 

c 


KRKC 11=1 
KKKC 2 I = 8 
P * 5. 

MR ITE< 6, 1000 1 

CALL PLUTM Y( XUOtfN , YACR OS , KKK * P J 

WRI TE! 6, 1U10 1 

RETURN 

1000 FORMAT! 2HPT, SOX, 24HBLADE SURFACE VELOCITIES* 

1010 FORMAT i 2HPL ,37X,63HVELuCl TV! «1 VS. MERIDIONAL STREAMLINE DISTANCE 


UMI DOWN 

THE PAGE /2HPL/ 




22HPL , 5CX, 

5 OH + 

- 

BLADE 

SURFACE 

1 . 

BASED 

ON 

32HPL , 5 CX * 

5 OH* 

- 

BLADE 

SURFACE 

1 * 

BASED 

ON 

42HPL, 5CX, 

5UHX 

- 

BLADE 

SURFACE 

2 , 

BASED 

ON 

52 HP L , 50X t 

SOHO 

* 

BLADE 

surface 

2 , 

BASED 

ON 

62HPL, 5U, 

5 OHS 

- 

BLADE 

SURFACE 

3* 

BASED 

ON 

72HPL t 5 CXt 

50h- 

- 

BLADE 

SURFACE 

3 t 

BASED 

ON 

B2HPL , 5GX, 

30HJ 

- 

BLADE 

SURFACE 

4 i 

BASED 

ON 

92HPL , SCXf 
END 

50H! 


BLADE 

SURFACE 

4, 

BASED 

ON 


MERIDIONAL COMPONENT / 
TANGENTIAL component/ 
MERIDIONAL COMPONENT / 
TANGENTIAL COMPONENT/ 
MERIDIONAL COMPONENT/ 
TANGENTIAL COMPONENT/ 
MERIDIONAL COMPONENT/ 
TANGENTIAL COMPONENT! 


SUBROUTINE VEL 

Va CALCULATES DENSITIES AND VELOCITIES FROM THE PRODUCT OF 
DENS ITT TiMtS VELOCITY 

COMMON SkW* E TER, tfcfO,L£*U2| ,MEfU2t 

COMMON /AUKRhO/ A (2000 , 4 1 » U! 20001 . K (2 0001 p KH UI2333 1 
COMMON /INP/ GAM ■ AR ■ T I P ,RH01 P , W T*"L t RTF LSF , OMEGA , QRF * BET A I , BET AO, 
LMB I , M BJ ,M B 12 T MBO 2 , MM , NbB I , m L ,NRSP , MR 1 50 1 t RMS P 1501 , BES PC 50), 

2BL OAT . A AN OK , ER SOR t StftFN, SLCRO . I N TVL # SURVL 
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CUMMJN /CALC UN/ MB (Ml ,MB iPl,M30Ml,MBQPl,Mbl2Ml .MBI2P1, MBD2H1, 
1M802P 1,MMM1, HM 1 > HM2 ,HM 3»HT ,D TLR ,OMLR i PI TCH ,C P , EX PON ,TWW , CPT IP, 
2TGR0G, raj , TBU, LAMBDA t TWL, I THIN, I TMAX ,N1P ,J NS (A) , BV (4) ,MV ( 100 ), 

31 VJ 101), I f VI IOC, 4), TV ( 100,4) .OTOMVtlOO ,4 I ,8ET AV ( LOCI ,4 i , 

4HH( 100,41 ,UTUMHI lu0,4l ,BETAH(100,4) ,RMHU00,4I ,BEH(10D,41, 

5KH( 100 1, BE( 1UU) ,DBUM( 100 ( , SAL (100) , AAA (100) 

CUHMUN /RHOS/ KHeiHB(i<J0,4) ,KH0V8 (100 .4 I 

OIMENS ION HWLRM1 100,4) ,WWCRTU00 ,4) .SURFLUOO ,4) 

U IRENS ION Wl 2000) ,RWMI 2000) .BE TA 12000) ,WMB (LOO ,4 ) , WT H ( 1J 3 , 4 1 , 
1XD0WNI Eim.YACROSI 800) 

EQUIVALENCE ( A ! I > 1 J , W{ 1) ) , ( A ( 1 ,2),RWMU)1 , ( A (1 ,3 I , BET A ( 1 ) ) , 

1< Al 1, 4 ), WMB( 1) I ,IA( 401,4) , WTb(l) ) , (A(801 ,4) ,XOONN(1I) , 

2! K( 1), VACR0S1 1 1 > 

C 

C V ELP CALCULATES ALONG VERTICAL HE SH LINES WHICH DO NUT 
C INTERSECT BLADES 

C 

ENTRY VtLP IF IR ST, LA ST, UPPER, LOWER) 

INTEGER BLOAT, AAN£3K,ER SDR, STRP N, SLCRU , 5UR V L , AA T£NP,SURF,Sl)HFBV , 
1FIRST, UPPER, UPPR8V,$1 .ST.SRW 

REAL X,KAR,L AHBDA , LMA X ,MH , HLE ,MR ,MSL , MSP , MV , MV I HI 
IF (FIRST .GT .LAST ) RETURN 

IF (FIRST. EQ.l. AND. (NTVL.GT.O) WRI TE (6,1 000 i 

IF (FIRST. EQ. 1) RELER = .0 

DU 20 JM=FUST,LAST 

IPU = IV( IM I 

1PL = IPU*NBSI-1 

TWLMR = 2.*OM£GA*LArtBOA-(OMEGA*RM(IM> ) 4*2 
LtR( 1 )=4 
DO 10 IP = IPU , IPL 
C UENSIY CALL NU. 4 

CALL DENSTY! H< IP 1 ,RHU< IP I , AN S , TWLMR ,CPT I P , EX PON, RNUI P » GAM, AR, T IP ) 
10 MI IP I - ANS 

IF I 1NTVL .LE .0) GO TU 20 

WRITE 16,1010) IH, IW(IP) ,B£TA(IP) ,IP=IPU,IPL) 

20 CONTINUE 
RETURN 

VELU0 CALCULATES ALONG VERTICAL HESH L I Nt S WHICH INTERSECT SLADES 


ENTRY VELBB l FIR ST, LA ST, UP PER .LOWER) 

IF ( F IKST .G( .LAST) RETURN 
IF ( FIRST .Nc .MB I ) GO TU 30 
SURH.I MSI, 1) = 0. 

SURFL(HBI,2) = 0. 

SURFL I MB I 2, 3 ) = 0. 

SURFL I MSI 2,4) = 0. 

30 DO 70 IM=FIHST,LAST 
I TV U = lTVI IM, UPPER) 

ITVL = ITV( EM, LOWER) 

1PUP1 = i P FI IM, I TVU) 

IPLM1 * IPF( IM, I TVL I 

TWLMR a 2 . *J ME6A*LAMS0 A* ( UMEG A* KM ( I M) 1 **2 
HCR = SURT ( TGKU&4TIP* ( 1 . - TWL MR/C PT1 PI ) 

IF I I TVL ,LT. ITVU) GO TO 50 
C ALONG THE LINE BETWEEN BLADES 
LERI 1 ) 3 5 

DO 40 iP = IPUPl, IP L Mi 
C DENSTY CALL NO. S 

CALL DENSTY! Wl IP l,RHO( IP I , AN S , TWLMR ,C P T I P , EXPON , RHCH P , GAM, AR, T IP) 
40 Wl IP) = ANS 

IF ( IN TVL .LE .0) GU TO 50 

WRITE (6,1010) IM, ( WCIPI ,BETA(IP) , IP*I PUP1 ,1 PLM1 ) 

C QN THE UPPER SURFACE 

50 RHUB = RHOVBI IM, UPPER ! 

LERI 1 J =6 

C DENSTY CALL NU . 6 

CALL OSySTYI WNB ( IM , UPPER ) .RHUVSIIM, UP PE RJ ,ANS , TWLMR ,CPT I P, 

IEXPON, RHU IP, GAff.AR ,TtP 1 
WMB( IM, UPPER ) = ANS 

WWCRMI IM, UPPER ) = WHB 1 I M, UPPER I / WC« 

IF ( IM .Eg .MB 1 .DR. 1 1 M. E U. MB 1 2 . AND .UPPER . EQ. 3 ) > GO TO 60 
DELTV = TV! JM-1,UPPER»-TV(IM, UPPER) 
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IF ( IM.ED.MdOP L.ANO. UPPER* EQ.3I Dc LTV = OE LTV- FI T CH 
5URFU 4H, UPPER I ** SUftFL!IM-l , UPPER* * SORT U MV l IM1 -MV t I M-l M **2 + 
II DEL TV ♦(Art! IMl*ftHI 1H-1U/2. 1**21 

60 RELER * AHAX11RELER,ABSI ( RHOB-RH DV0II M f UP PERI J / RHOV BI I M f UPPER! II 
C ON THE LOwfcR SURFACE 

ft HO 6 * RHGV&l IH,LGW£fcl 
LER! 1 1 =7 

C OEMST* CALL NO* 7 

CALL DENSTYI WM0UM,LOWtftl ,RHOV8! lM,LOW£ftl , ANS , TWLHR , CPT I P, 

1EXPUN • ft HO IP t GAN t AR t TIP I 
WHB i IHiLOWEft ) = ANS 

W WCRM ( 1M, LOWER I « WMii! iMjLOWERI/lCR 

IF I I H *E3 *M0 I *0 R * ! I M, E Q . MB 1 2 . A ND * LOWER* E 0- 41 ) GO TO 70 

DEL TV = T V 1 1 M* l , LO HER 1 - T V 1 1 M r LOWE R I 

SURFLt IM* LOWER) = SURF LU M-l , LOWER} * SQRT U HV U K) -HV 1 1 M-l >t **2 + 

li OELT V*! RMI I HI +RH( IM-1 1) /2. 1**2 1 

TO ft EL Eft = AM AX 14 RELER *AB SI ( RHGB—RHOV B1 1 H , LOWER! I / RHOV Bt I H, LOWER) I } 
RETURN 

VELSUft CALCULATES ALONG A B LAOL SURFACE 

ENTRY VELSUft 
DO 90 SUftF-1,4 
IMSS - IM ST SURF I 
IF! iH5S.E0.0i GU TO 90 
DO 80 IHS-1, IMS5 

TWLHR - Z **0 MEGA* LAMBDA** 10MEG A*RNH 4 IHS ,SURFI I **2 
WCR = sqrtitgrog*tip* {I.-TWLHk/CPTlPH 

RHUB “ RHOHtil IHSt SURF » 

LER* 1MB 

OENSTY CALL NO. 0 

CALL OENSTYi WTB ( IBS, SURF ) *RHOH B 4 IH S, SURF * ,ANS ,TWLHft,CPUP, 
IEXP0N*RB3 IP , GAM * AR * T LP 1 
WTB4 IHS, SURF * « ANS 

WWCRT! IH$,SURF1 - WTBf IH $ , SURF 1 / wC R 
00 RELER * A HA X 1! ft EL E ft , A & S ! 1 RB U0- ftHOH B (IHS, SURF II/ RBOBB TINS, SURF) J J 
90 CON 1 IN UE 

IF (RELER.LT ..0C1I IE NO * IEND+1 
WRITE! 6, 10801 ITER, RELER 

WRITE ALL BLADE SURFACE VELOCITIES 

IF 4 SURVL *L£ *01 RE TURN 
WRITE! 6, 1020 1 

WRITE! £* 10 AO J mVUm i WH0! IH.l I ,B£ TAVUMtll ,5URFL! I M , L 1 , 
iWWCRHt IN, 11, WMBUH,2) ,aETAVUH,2 I , SURF L UN, 21 ,WWCRMUM,2I, 
2IM=MBI,M&0 1 
WR Utl t, 10 JO J 

WR 1 TE( i f 10 AO 1 (MV< I Hi ,WMB( !M,31 , 3£TAV(IM,3 1 ,$URFL( IM,3 i , 

1WWCKHI lM r 31, WM&( IH,41,8ETA V(IM,4l , SURFLUH.9J , HWtRH t I M , 4 * , 
2IH~MBU,MBU2 1 
WR ITE! tr 10 501 
DU 1UU SURF*!, 4 
IHSS * IMS! SURF! 

IF! IHSS.Ll.Li GO TO LOO 
WR I TEC 6, 1060 I SURF 

WRITE! i t 10701 (MHI IHS, SURF 1 , WTB ( IH S , SURF 1 *BETAN! IHS ,SURFi , 

1WWCRTI IHS, SURF 1, IHS^l, IHSS1 
100 CONTINUE 
RETURN 

100 0 FORMAT UH1////40X,34HVE LUC I TIE S AT INTERIUR HESH POINTS// 1 
10 10 FORMAT ! IHL , 3HlH & , 1 3, 5! 24H VELOCITY ANGLE (DEG) }/ 

H SX,5( G15.4, F9.21 1 1 

1020 FORMAT! 1H1////1GX, IH* , 10X,49H SURF ACE VELOCITIES BA $ ED ON MERIDIONA 
1L COMPONENTS, 43X, lH*/16X,tH* ,53X,tB* ,5 6X,1 H*/16X ,1H* , 19X , 15HBL ADE 
2SURFAC E l f HX, IH* ,2OX,i5HBLA0e SURFACE Z *21 X, t H*/7X ,1HM,BX, IH*, 2! 3 
3X , 0HV £LQC i TY , 3X , 2 3HANGLE (DEG J S URF . LE NOTH ,5 X,5 HU/ WC Rt 6X , 1 H*» 3X H 
1030 FURMAH///1EX, IH* , 19X, ISHBLAOE SURFACE 3 ,19X,IH*,2 0X,15HBLADE SURF 
1ACE 4, 21X, lH+/7X,lHHiBX,lH* t 2 13 X , B HVE LOG I T Y, 3X ,23HANGL El DEGI SJftF, 
2 LENGTH, 5X,SHW/WCft,6X T IH* ,3XH 

1040 FORMAT (IH ,013. 4, 3H * ,G 12 . 4 ,F 9. 2 ,2G1 5.4 ,6B * ,GI2.4,F 9.2, 

L2G1 5* 4 , 3H *J 

1050 FORMAT 1 IH 1// //3X * 49H SURFACE VELOCITIES BASED ON TANGENTIAL COM PONE 
1NT S } 

1060 FORMAT I //22x f lSHBLADE SURFACE , 1 1 /7X t \HH UDX , 0HVE LOC ITY, 3X , 10 HANG 
1L El DEG }, 3Xt 5HW/WCR 1 
1070 FORMAT! IH , 2G13. 4,F 9, 2,G 15* 41 

1080 FORMAT U4HL ITERATION NO. , 1 3 • 3X ,36H MAXI MUM RELATIVE CHANGE IN DENS I 
ITY ", G il * 4 1 
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SUBROUTINE GLCO 

BLCD CALCULATES BLADE THE TA COORDINATE AS A FUNCTION OF M 
CuMMGn SKW. { I Eft . I END ,LfcR(2) ,NER(21 

CO MM Jo / 1 nlP / GAM.AK.TI P.RHUI P.WTFL.kTFLSP, OMEGA. ORf, BET Al , BETAS, 
IM0I,M0U,N012,MBU2,MM,N801 , NO L . NRSP , MR 1 50 1 ,KMSP (50) , BESPtSD 1 ,• 

2BL DAT , AAMOiC , £R SOK , STRF « , SLL RD , I N T V L , SU R VL 
COMMON /CALCUN/ MOIMl.MBIPl ,M0OM1 , MBUP l , MGI 2 Ml ,MBI2Pl , H0U2H1, 
1MB02P l.MHMl, HM1,HM2.HM3,HT .QTLR.DMLk ,PI TCH »C P , EX PUN ,T WW , CPT fP, 
2TGKU6, TBl, T»U, LAMBDA, I WL.l TMI N,( TMAX ,(Mt P ,1 MS I A 1 , BV 14 1 , MV ( 100 ) , 
31V ( 101 ), IT VI 1U0.4) ,TVU00,41 ,0 TO MV 1100 ,41 > BET AV ( 13 0 ,41 , 

4MHI 100, 41, DTuMHI 100,41 ,6E TAN ( 1 00 ,4 1 ,RMHU0O,4l , 6EH 1 10 3 , 4 I . 

5ftH( 100 l,BE 1 100 1 , OB DM l 100 1 , SA LI 1 001 ,AAA(1001 
COMMON /GEUM IN/ CHORD I 41 *ST6R(4) ,MLE (41 ,THLE (4 1 ,RM1 (4) ,KM0<4), 
IKK 41 • RlM 4) , BE T I < 4) ,gE 10 1 41 ,NSPI (41 , MSP (50 ,4) , THSP (50 ,41 
COMMON /ULCOCM/ EMI 50,41 , INI T (41 
ENTRY tit MM , THETA »DTOM, INF 1 

INTEGER dL DA T , AANDK »ER SOR , STKF N i SLC KD , SDK V L , A A T E HP , S U RF , SU KF BV , 
IF IK ST, UPPER, UP PR BV, SI, ST, SR W 

REAL X , KAK ,L AMBOA ,LMAX ,MH , MLE , MR ,MSL , M SP , MV , MV 1 Ml 
REAL M.MMLE.MSPMM.MMMSP 
SURF= 1 
SIGN:, l. 

GU TO 10 

ENTRY BL 2(M » THE TA »UT0M« INF I 
SURF* 2 
SIGN* -1. 

60 TU 10 

ENTRY BL 3IM , THE TA ,0 TOM, 1 NF 1 
SURF* 2 
SIGN* 1. 

GU TU 10 

ENTRY a. 4IM , THE TA ,0 TDM, INF ) 

SURF* 4 
SIGN* -1. 

10 INF* 0 

NSP* N SP 1 1 SURF 1 

IF I IN 1TI SURF). EG. 131 GU TO 30 
1NITISLRFI* 13 

INITIAL CALCULATION UF FIRST AND LAST SPLINE POINTS ON BLADE 

AA = 6 £T 1 1 SURF 1/57,245779 
AA * S INI AA I 

MSPII.SURF) * R 1 1 SURF )*I1«-S1GN»AA1 
OB = SORT! l.-AA**2) 

THSP< 1, SURF 1 = S1GN*00*M I SURF I /RHMSURF1 
BETH SURF! * AA/BB/RMM SURF 1 
AA = G ETUI SURF 1/57. 295779 
AA * S IN I AA 1 

MSPlN SP, SURF I * CHORD! SURF l-KDlSURfl *11. *SI6N*AA1 
00 * SORT I U-AA**Z) 

THSPI N SP , SUR FI = STGRI SURF I * SI GN*BB*RO ( SURF I /RMUISURF ) 

BETUISLRF) * AA/0O/RMO l SURF 1 
00 20 IA = 1,NSP 

MSP I I A, SURF 1 » MSPIIA.SURFKMLE (SURF) 

20 THSPI I A, SURF I* ThSPU A , SURF 1 + THLE I SURF 1 

CALL SPL 4 221 MSP! I, SURF 1 , THSPI I ,SURF I ,BETf (SURF I , BE T 0( S URF 1 , NS P , 

I AAA, EMI 1, SURF I I 
IF(OLOAT.LE.O) GO TO 30 
IF (SUKF.EU.il WKI TE 16,1000) 

WRITE! E, 1010 I SURF 

WRITE (6, 1020) (MSP (I A, SURF I , THSPI I A ,SURF I , A A A ( ( A1 ,EM( I A, SURF) , 
1IA=1,NSPI 

0LAUE CUUR01NATE CALCULATION 


30 KX * 2 

IF IM.CT.MSPI 1, SURF I) GU TO 50 
C 

C AT LEADING EUGE RADIUS 
C 

HMLE* M-ML6ISURFI 
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if i MMLc.LT *— dnlr ) GO TO 90 
MML6* AHAXII0*,MML£* 

TH£TA= S0HT{ MKLE+l 2»*ft I f SURF J ^MMLE 1 1 #51 GN 
IF I T HETA *E0 #0* I GO TO 40 
RMM = RHSURFI-MMLE 
OT DM* RMM/THETA/KMUSURFI 
THETA* THE TA /KM 1 1 SURF J +THL£ { SURF I 
RETURN 
40 INF* I 

DTOM = 1 *E 1Q+5 IGN 
THETAs THLE15URFI 
RETURN 

ALONG SPLINE CURVE 

50 If t M .LE *M SP IKK .» SuftFI I GO TO 60 
IF ( KK *GE *N$P I GO TO TO 
KK * KIU1 
GO TO 50 

60 S* MSP IKK r SURF HMSPlKK-1 , SURF I 
LMKMl* EM C KK** 1 * SURF ) 

£HK* EMlKK t SURFI 
MS EMM* M SP IKK, SURF *-M 
MMH SP - M~M$P(KK~i,SURFl 
THK* T TSP 1 KK » SURF I /S 
THKM1* THSPtKK- 1,SURF 1/5 

THETA* £MKM JL+MSPMM** 3/6* /5 + E MK+MMNSP++3/6* / S * I T HK-EMK+S/6* l* 
1MMM SP + ! THKH1-EMKM1+S/6* 1 +MSPMM 

OTOM* -FMKM1+MSPMM++2/2./S * E MK+HMMSP++2/2. /$ * FH K**T HKM1 - I 6MK- 
1EMKM 1 *#S/ 6* 

RETURN 

AT TRAILING EDGE RADIUS 

TO CMM* CKJRD( SURFHMLE* SUKFI-M 
IF !CMM*iT*-UMLRI GU JQ 90 
CMM a A MAX 1(0*, CMM I 

THETA* 53R f{ CMM* ( 2 *+R D( SURF * -C MM* J + S1 G N 

IF ( THETA *£Q *0* I GO TU 60 

RMM- RQ< SURF 1—CMM 

DTUM = -RMM/ THETA /RHQI SURF } 

THETA = STGRISURF THE TA/RM0( SURF* ♦THLEf SURF I 
RETURN 
SO INF- 1 

Of DM - - 1 *E10* SIGN 

THETA* THL£i SURF I + STGR I SURF * 

RETURN 

ERROR RETURN 

90 WRITE! £, 10301 LER 121 , M , SURF 
SfOP 

1000 FORMAT ( INI, 1 3X , 33H0LA ilE DATA AT INPUT SPLINE POINTS* 

1010 FORMAT I 1HL , 1 ?X , 16HGLA0E SURFACE,! 41 

1020 FORMAT | TX t 1HM , iOX, 5H THE TA ,1 OX * lOHOERI VAT I VE , 5X ,10H2ND DERIV* / 

H 4GIS- sn 

1030 FORMAT U4HLOLLO LALL NO. ,I3/33H M COORDINATE IS NOT WITHIN BLADE/ 
14H H *,G1 4*6, 10X, GH SURF *,G14. 6) 

END 


FUNCT ION IPFUMpI T* 

COMMON /IHP/ G AM , AR , T I P , RHOt P , H TF L * WFF L$P T OMEG A , UR F , SET A l , BET AO * 
IMG 1 1 M 00 p i 0 1 2 , MOD 2 p MM , N8 8 1 , N0L , NR SP , MR ( SO* ,RM$P(5€H ,flESP(50 I, 

26L DAT , AANUK, £R SOR, STRF N , SLCRO , 1 N TVL , SURVL 
COMMON /CALC ON/ MB IM 1 , MS IP i , H8QMI ,MQ 0P1 , MB l£ Ml , Mil 2 PI , MS02M1 , 
1MB0 2P 1 jMIM 1 , HM I T HMZp HM3,H T ,0 TLR ,DMLR *PI TCH ,C P > EXPOM ,TWW , CPT I P* 

2T GROG p TbI f TGOtLAMSDApTWLplTMINp I TMAX p NI P ,1 MS 1 4 * , 0V (4 } , MV 1 1(JQ I , 
3tVi IQI ), IT Vi 100,4* ,FV( 100,4* ,DTOMVUOO ,4i ,0ETAVmO,4J, 

4MH( 100 , 4* , DTDMH 1 TOO, 4 * ,RETAHU00,4l t RMH(100,4l ,B£H( 103 ,4* t 
5RM( 100 ** 0£| 1001 fDOUMllOO* ,SAL!1GQJ ,AAAUOO* 

IP F = IV( IM* + fT- ITVIIM,! *- I TVU M,3|-10000 
IF! irt.LT.MBI 2*0R*IM*GT*MBUI RE TURN 
IP F - IV! 1M> +IT~ITVlirt,U 

ifi 1T*6E. ITVUM.3I I IPF * I PF - 1 TVU rt ,3 1 * 1 T VII H t 4 * + 1 

RETURN 

END 


Subroutine MHORIZ 


Subroutine MHORIZ calculates the m -coordinates of intersections of all horizontal 
mesh lines with a blade surface. It locates points of intersection, by checking the ITV 
array (see main dictionary). E ITV changes between adjacent vertical lines, there must 
be a horizontal mesh line intersection between those vertical lines. ROOT is called to 
calculate the m -coordinate of the intersection. The input arguments for MHORIZ are as 
follows: 

MV array of m -coordinates of vertical mesh lines 

ITV same as ITV of main dictionary, but for a particular surface 

BL subroutine giving blade 6 -coordinate as function of m (BL may be BL1, BL2, 
BL3, or BL4 in the calling statement. These are the entry points of BLCD. ) 

MBI value of EM at first vertical mesh line to be checked 

MBO value of EM at last vertical mesh line to be checked 

ITO value of IT at the origin of coordinates at leading edge of front blade 

HT mesh spacing in 9 -direction 

DTLR tolerance in 6 -direction 

RODE code variable to indicate whether blade surface is upper or lower; RODE = 0 for 
upper blade surface, RODE = 1 for lower blade surface 

MRTS integer switch indicating infinite slopes at leading or trailing edge of a blade sur- 
face 

The output arguments for MHORIZ are as follows: 

J counter indicating current number of intersections of horizontal mesh lines 

with a given blade surface 

MH m- coordinate of intersection of horizontal mesh line with blade 

DTDMH slope d 6 /dm where horizontal mesh line meets blade 

The internal variables for MHORIZ are as follows: 

IM vertical mesh line number 

ITIND counter of horizontal mesh lines which intersect blades between two consecu- 
tive vertical mesh lines 

MVIM MV at left end of horizontal interval 

TI 0 -coordinate of horizontal mesh line which intersects blade 
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SUBKuyi 1NE MHURIZ1HV, l TV,SL ,HBI ,HBOpl TOiHT ,UT L ft , KOBE , J , MH, 01 DMH, 
lMRTS I 
C 

C HHQRIZ CALCULATES M CQQRUINArES OF INTERSECTI QMS OF ALL HORIZONTAL 
C MESH LINES rfiTH A BLAOE SURFACE 
C RODE = 0 FOR UPPER 0LAUL SURFACE 
C RQQE = 1 FOR LOrtLR BLADE SURFACE 
C 

CUMHUfM SR W , l TE R , I END , LLR 1 2 J f Nt R ( 2 } 

DIMENS iOH MVI 100* * I TV! 1UQ* T MH ( 1 0 0) «0 TOW U (ill 

INTEGER 6LDAT,AAN0R,tR SUR, STRFN, SLCRO r SURVL,AATEMP ,S J RF , SU RF BV , 
IFIRSTt UPPER* UPPRBVt Si , ST# SKK 
REAL K ,KAK il AMUOA , LMA X , MH , ML5 r MR ,M5L , MSP * MV i MV I HI 
REAL H V 1M 
EXTERNAL liL 

If tMaip.GE.MBO I RETURN 
1H= MB I 
10 IT I N Q = 0 

ZO IF f IT VI IM + i J- I T VI IMI-I UNO* 30,40,50 
30 J= J*1 

Tl~ FLOAT l IT Vi 1M + U-1 TU- I TI ND + KOBE I *BT 
IT |ND = 1T1ND-I 
MV IM = MVI 1M I 

IF IMRTS.EU.l) MVIM = MVI M+ ( MV U M* 1 J -MVI Ml / 1 33 0* 

CALL ROOT tM ViM,MVtlM+ U til p 3L rDTLR, MHI JI fOTOMHUH 
GO TU ZO 
40 IM~ IM+1 
MR T S =5 0 

IF I IM .kg *M60 1 RETURN 
GO TO 10 
50 J- J + l 

T I— FLOAT < I T Vi IM|- ITU+I TI NU + kGOE J*HT 
IT IND = IT INO + 1 
MV IM = MVI 1M i 

IF tMRIS.bg, H MVIM = MVI M+ { MV ( l M* L ) —MVI MJ /I OD 0 + 

CALL ROOT IM VIM v MVt £M^1) r II ,BL fOTLK* MH ( J) pDTOMHt J) I 

GO TU ZO 

END 


Subroutine DENSTY 

Subroutine DENSTY calculates the subsonic relative velocity W and corresponding 
density p that result in a given value of the mass flow parameter pW. This is done by 
using equations (B5) and (B6), which are an algorithm based on Newton's method. 

E the value of pW is too large, there is no solution. In this case an error message 
is printed, W cr and the corresponding density are calculated as output, and the program 
continues. Thus, it is possible to get an approximate solution even though there may be 
one or two points with too large a value for pW. The input arguments for DENSTY are 
as follows; 

RHOW pW 

RHO initial estimate for p (pj n may be used) 

TWLMR 2wX - (o>r) 2 

CPTIP 2c T! 

P m 

EXPON l/(y - 1) 
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RHOIP 

GAM y 

AM R 

TIP 

VTOL convergence tolerance on relative change in W 

The output arguments for DENSTY are as follows: 

RHO p 

VEL W 


The internal variables for DENSTY are as follows: 

RHOT newly calculated estimate for p 

RHOWP d(pW)/dp {eq. (B4)) 

TEMP ( T / T ^)(2-r)/{?-l) 

TGROG 2yR/(y + 1) 

TTIP T/T! n 

VELNEW newly calculated estimate for W 


SUBROUTINE OENSTriftHtJW^RHQrVEL, TWlMRfCPTl P , EXPCiN, RHtJI P , GAM, Art, T IP) 
C 

C DENSTV CALCULATES UENSITV AND VELOCITY FR CM THE WEIGHT FLOW P ARAMET Eft 
€ DENS 1 TV TIMES VELOCITY 
C 

COMMON SRWiITeRafeNOiLeRiZf *N6RI2* 

VEL =* ftH3 W/ft HU 
If lV6L.Nt*0.> GO TO 10 
RHO = RHO TP 
RETURN 

10 TTIP = i.-t VEi**2+TWLMftl /CPTIP 
tfi TT tP.LT*0 + j GO TO 50 
TEMP * TT IP* *4 EXPliN-l* i 
RHOT ~ ftHOlP*TEHP*TTIP 

RHOWP * “ VEL** 2/GAH*RHGl P/Aft+TE HP/T I P*- RHOT 
£Fi RHQ mP *LE * 0* J GO TO 30 
VEL NEW * VEL +T ftHUW-RHQ T* VEL I /RHOWP 
IFt ABSlVELNtW-VELI/VELNE W,LT- ,000i» GO TO 1 0 
VEL = VELNEW 
GO TO 10 
20 VEL = VELNEW 
ft HQ = RHUW/VEL 
RETURN 

30 TGROG - 2**&AM*Aft/£GAH+U I 

V£L * SORT 4 T GRO&*TlP* 1 1 * - T WLHft/C PT I PI 1 
R HO * RHQ IP* 41 *-{ VEL**2*TWLMR)/CPTJPJ**EXPON 
RWMOftW * ft HU W/ftHO/ VEL 
NEfti U - NERI1J + 1 

WRITE! e, 1000 i LEfttU.NERUt ,ftWHURW 

lF{NEft(iT*E0*50» STOP 

RETURN 

1000 FORMA T £ 1 6HLU Ett ST V CALL MU, ,E3/9H NEftQI *,I3/1QH RHQ*W IS f P7* + * 
134H TIMES THE MAXIMUM VALUE FOR RMO*W) 

END 
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Subroutine ROOT 


Subroutine ROOT finds a root for f(x) = y by Newton's method. The function f(x) 
must be defined on a specific interval [a,b] . The values of f(x) are calculated by 

another subroutine (FUNCT). 

k+1 k 

The value x ' is determined from x by 

x k+ i = y - fo k ) + x k 

f’(x b ) 

The first value of x is x® = If x b+ * is not in the interval, if f'(x b ) = 0, or if 
f f (x k ) = the interval [a, b] is scanned to see if a suitable starting value of x for 
Newton's method can be found. If a root cannot be found in 1000 iterations, a message 
is printed, and the calculations are stopped. 

Subroutine ROOT requires that f{x) be calculated by a FORTRAN subroutine sub- 
program {FUNCT). Any name may be chosen for this subroutine. In TANDEM, FUNCT 
is either BL1, BL2, BL3, or BL4. The calling sequence is 

FUNCT (X, FX, DFX, INF) 

These arguments are defined as follows: 

FX f(x) 

DFX f’(x) 

INF used to indicate an infinite derivative: 

0 if f'{x) is finite 

1 if f'(x) is infinite 

The input arguments for ROOT are as follows: 

A a 

B b 

Y y 

FUNCT external subroutine to calculate f{x) 

TOLERY tolerance on solution (x is accepted as a root if jf(x) - yj < TOLERY.) 
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The output arguments for ROOT are as follows: 


X value of x such that f(x) = y 
DFX f'(x) 


SUBROUTINE ROOT* A t fl r V , F UNO T , TOL ERY t X , DF X 1 
C 

C RCGT FINOS A ROOT FOR IFUNCT MINUS Y) IN THE INTERVAL f A T B 1 
G 

COMMON S R W t ITER, I END, LER f 2 ) , NER t 2 J 
INTEGER SRW 

IF (SftW.EQ.2U WRITFI&, 1000) A * B , Y , TOLER V 
TGLERX= tB-A)AQOO. 

A82 = t A + R ) /2* 
i = 0 
X = A 

10 CALL FUNCTUtFXfOFXtlNF) 

IF fSRW.EQ.2l> WRITE(6, 1010) I f X f F X * DF X , 1 NF 
IF I ABM Y-FX) .LT .TOlERY ) RETURN 
IF l WOE, 1000) GC TO 30 
1 = 1 + 1 

IF t INF *NE .0 .OR * DFX-EQ.Q.) GO TO 20 
X= tY-FXl/DFX+X 

IF (X.GE.A .AND- X-LE.B) GO TO 10 
X = A+TGLERX*FLQ Alt l ) 

!F( f*EQ. 1 ) X = B 
GO TC LO 

20 IF iX.LT.AB2S X=X+TOLERX 
IF (X.GE.AB2J X=X-TOLERX 
GO TO 10 

30 WRETEtb, 1020) LERt2),A ( 8,Y 
STOP 

1000 FORMAT ( 32H 1 INPUT ARGUMENTS FDR ROOT — A =GI 3. 5 , 3X , 3HR >,GI3.5, 

1 3X , 3HY = t G13.5,3X,BHTOLERY *,GI3.5/17H ITER. NO- X,!7X, 

22HFX,iBX r 3HCFX t 10X,3HtNF) 

L010 FORMAT ISX, I 3 f G1 6 . 5 1 2Gi B* 5 , I 6 ) 

1020 FORMAT (14HLR0GT CALL NG»,13/47H ROOT HAS FAILED TC CONVERGE IN 10 
LOO ITERATIONS/*** A -,G14.6 f I OX , 3HB = * G 14 . 6 , 10X t 3HY =,G14.6) 

END 


Subroutine SPLINE 


This subroutine is based on the cubic spline curve. SPLINE solves a tridiagonal 
matrix equation given in reference 9 to obtain the coefficients for the piecewise cubic 
polynomial function giving the spline fit curve. SPLINE is based on the end condition 
that the second derivative at either end point is one -half that at the next spline point. 
The input variables for SPLINE are as follows: 


X array of ordinates 

Y array of function values corresponding to X 

N number of X and Y values given 
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The output variables for SPLINE are as follows: 


SLOPE array of first derivatives 
EM array of second derivatives 


If Q = 13 in COMMON, input and output data for SPLINE are printed. This is useful 
in debugging* 


SUBROUTINE SPLINE [ A , Y , N t SLOPE t E M) 

C 

G SPLINE CALCULATES FIRST A sMU SECOND DL R i V A TI VE 5 AT SPLINE POINTS 
C END CuNDITIJN - SECOND DERIVATIVES ARE THE 5 A HE AT END POINT AND 
C ACJACENT PUINT 
C 

CUMMuN W/ttOX/S( 1001 , A < LOO I »B I100J ,C 1 1001 ,F (1 OQI r Wt 10J ) rSB( 100 \ , 

1GI 200 } 

DIMENSION XI NJ t VIM ,£H t NJ , SLUPfc ( N| 

INTEGER 0 
DO LG ! = 2 t N 
10 St I >=Xt i I — XI. 1-1 * 

ND=N- i 

IH .MU ,LT *21 GO TU 30 
00 20 1=2 ,NJ 

At I ) = 5 t I ) / 6 . 

fctl n = t Si 1 J+St i + IM/3. 

Ct I J*st 1 + H/6, 

20 FI i * = t VI l+l|-Yt 1 H /sii+il-i Ytn-YU -m/st H 

30 At N \ = - , 5 

at i J = 1 * 

Bt H )=1 . 

Ct 1 I * - *5 
Ft 11 = 0* 

Ft N 1 = 0* 

H{ L 1 = 0(1) 

sat i i^ct li/^ii 

Gt 11=0 * 

DU 40 1 = 2, N 

mi ihiiiii-Ai u*sat i - u 

sti t n = ct n/wtn 

40 G( n = [ Ft 1 l-A t I >*GC i-i I J /*U i 
EM I i» J = Gi M J 
DU ^0 I = l , (V 
K^N + I- I 

50 EMt&| = GliU-SBlKl*fcM(K+iI 

SLUPtt U=-St 21 /6** t2.*L-Mt 1 J+EMI2 U * t Y12) -Y <1 I I /St 2 I 
DU 60 1 = 2 

6 0 SLuPEt IJ = Sl I l/6.*t 2.*EMtl I +E Kt I — 1 » I + ( Y < I ) - Y f 1 -1 J I/S ( II 

IF 4u.Eu.13J HR iTt t6 f LOGO) N • t Xt 1 1 , Y (II , S LOPE C II , E M t 11 , 1= I , N I 
RETURN 

100 0 FUkMAT tiX,L5HN0* DF Pul NTS = 1 1 3/1 0 X ,1 H X , 1 9X , 1 H Y » 1 9 X , 5 HS LD PE » 1 5X , 
1 2HEM1 / 1 4F20.Q I I 
END 
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Subroutine SPLN22 


This subroutine is the same as SPLINE, except that, for the end conditions, the 
slopes are specified. The input variables for SPLN22 are as follows: 

X array of ordinates 

Y array of function values 

YIP slope at first point 

YNP slope at last point 

N number of X and Y values given 

The output variables for SPLN22 are as follows: 

SLOPE array of first derivatives 

EM array of second derivatives 

If Q = 18 in COMMON, input and output data for SPLN22 are printed. This is use- 
ful in debugging. 


SUBROOTlut SPLN 22 IX, V, YIP, YNP.N, SLOPE ,EM 

c 

fc SPLN 22 CALCULATES FlftSl AND SECOND uEKl VAT! VE 5 AT SPLINE POINTS 

c end condition - derivatives specified at end points 
c 

COMMON U/BJX/Sl 100 » ,A 1 1001 ,B tlOO) ,C UOO) ,F 1100 ) ,W 1100 ) ,SB! 100 t , 

nil 2ou ) 

LHMfcNS lUN Ai N ) r (iN * t SLOPE 4 M 

mtGEK J 
00 10 I^rN 

io si i { n-x* i- 1 ) 

UU^N-l 

IH faG *LT *21 Gu 10 30 
DO 20 I=2*NJ 
A{ i/6* 

i j~t Si i i+si i + iR/3* 

Ci M*St I + i i/6* 

20 Ft i i=t vt i + d-yu M/ siH-iJ-ivin-vu-m/sui 

30 Ai H > = SIN)/ 6 * 

0< 1 l*Si 2J / 3. 

0( N } = St N )/ 3 * 

Cl ii = Si2J/6* 

H 1 1 = 1 Yi 2 l-Y ( ill /Si 2i- YIP 
Rn) = YW-( YiN >-Y(N-ii J /SIN) 

Wi 1) = 6UI 

SfcU H=Ci i 1 /*U l) 

G( U^Ft 1 I / W I li 

UO 40 1*2, N 

W< I l-d ( 4 i-Al i i * $b 4 1 ~ 1 ) 

$B 4 n=CiI)/rtUI 
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40 b( I 1=1 K i i-A( i )*G4 1-1 n /Mi L J 
tMiNI = GUM J 
00 b'O !=*2*N 
K-N + l- L 

5G tHl IU = G(K » — Sti 4 K 

SL^Pt i 2>/t>**I2.*EM(lUfc:Ml2n + ( V ( 2 > - V 11 H / S <2 \ 

00 6U 1^2 $ N 

6€ SLUPtl ll»S{ I l/6.M2.*ErtU |ttmi-inMYUi-Y{Wl l/SM } 

IF UU£0* LfcM *tK I f E 1 6 T 1 DO 0 ) N , I XU ) » V l ll t$ LOPE IM , E M 1) t 1= 1 * N I 
KET UKrt 

1000 FORMAT f2X t i5HNO* OF JOINTS = , I 3 /I OX ,1 HX , L SX , IHY ti 9 X ,5 HS LOPE » 1 5X t 
I2H£H/I 4F2U,B } I 
END 


Subroutine SPLINT 


This subroutine is based on the cubic spline curve, with the same end conditions as 
SPLINE. The cubic spline curve is then used for interpolation. The input variables for 
SPLINT are as follows: 

X array of spline point ordinates 

Y array of function values at spline points 

N number of X and Y values given 

Z array of ordinates at which interpolated function values are desired 

MAX number of Z values given 

The output variable for SPLINT is as follows: 

YINT array of interpolated function values 

If Q = 16 in COMMON, or if some element of the z array fails outside of the inter- 
val for the x array, input and output data for SPLINT are printed. This is useful in de- 
bugging. 
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SUBROUTINE SPLINT I X , ¥ ,N ,Z f MAX # YI NT ,D YU XJ 


C 

C SPLINT CALCULATES INTERPOLATED POINTS AND DERIVATIVES 
C FOR A SPL INE CURVE 

C END CONDITION - SECOND DERIVATIVES ARt THE SAME AT END POINT ANO 
C ADJACENT POINT 
C 

COMMON O/bUX/SI LOOl ,AUOui ( CI100I ,F Cl 00 1 103* ,$0C 10DJ, 

16C 100)»EM< 10 01 

DIHENSIUM XINl,r(Nl ,Z(MAX* , Yl NT! MAXI ,OYOXCMAX| 

INTEGER 0 

IP 4 MAX mLt .0 J RETURN 
ill = 0 

00 10 1 = 2 , N 

io si n=xc n-xi i- ii 

nu=n-* 1 

1PCNU .LT.2* GU TO BO 
DO 20 I - 2 , NO 
At I | = St n/6*c 

bi n = t si i *+si i + ii 1/3*0 

Ct I 1 = S 4 I + H/6.G 

20 PI i 1 = 1 YU + lirVU M/SU + 1IMYU I/S II I 

30 At N ) = — . 5 

01 11 = 1*0 
BIN 1 = 1*0 

C( 11 - ” . 6 

PI 11=0*0 

FIN 1 = 0.0 

Wl 11 = 0111 

SBC U=U ll/wll) 

GC u = o*o 
DU 40 l = 2,N 

rtl 1 l = BC 1 1“ A 1 I **StH I-ll 
sbi ii =a n/KU i 

40 G( 11 = 1 PI I l-AC i 1*GI l-l I I /rftl 1 
EM I N 1 = G( N J 
DU SO j = Z,N 
K*N+i- I 

50 EMIK1 = GIK *-$BCK»*EMCK+ll 
DU 140 [= 1 * M A X 
K-2 

IPU4 I HXi II I ?0t60»90 
60 Y INTI I 1- Y I II 
GU Tu 130 

TO IRZC U,bE*( 1*1*XC 1J-,1*X<2H l GO TO 120 
WRITE 16, 100 Q 1 ZU \ 

U = 16 
GU TU 120 
00 K =N 

IFIZC ll.LE.I l.l*XlNI-.l*X(JWm G U TU 120 
Wk ITt it ?, 10001 2111 
0 = 16 
GU TU 120 

90 IFIZlIHiMI 120,100*110 
100 Y IN T( i *=YlM 
GO TU 130 
UG K -K + 1 

fflK-N 1 9 C, 90 , eo 

120 YINKI) = bM(X-ll*CX(iO-2(I 11 **3 /6. / S I K| + EMI Kl CZ ( I I -X ( K-l* I ** 3/ 6. 
l/SIM + I YIK I/SIM-EMIRI *SU*/6. 1 * U< 1 1-XtK-ll 1 + t Y I K-l 1 / S I K I -EMI X- II 
2*SU 1/ t* )*t XIK 1-24 1 1 \ 

130 DYDXC Il=-ENt K-1|*C XiKl-Zl 111 **2 /2. 0 / S C KI + E MC KJ * ( Xf K-L 1 -Zt I 1 1**2/ 2. 

10/SIK I H YCR I -Y IX,- ill ZStKl-Ct MUU -t MtK-ll 1*5110 76.3 
140 CONTINUE 

MX A = M A X 0 1 N » M A X 1 

IPI U *Eg .16 1 WH I TEt b, 10101 N, MAX, IX ill , YUl ,21 II , Yl NT i I 1 , UY OX ( I) * 
11= 1, M X A I 
U - III 
RETURN 

1000 FORMAT I 54H SPLINT USE U PUR E XTRAPULAT f ON. EXTRAPOLATED VALUE = * 
1G14.6I 

1010 FORMAT i2X r 2lHNU. UP POINTS GIVEN = ,I3,30N, NO. CP INTERPOLATED PO 
1ISTS I3/10X, 1MX, 19X, lr*Y,16X, II HX- INTERPOL. ,9 X, 11 HY-I NT ERPOi*, 
2BX, 14H CYUX- I N T LRPOL . / i 6E20-BM 
ENU 



Lewis Library Subroutine TIME! 


This subroutine is part of the Lewis Systems Library. TIME1 gives the time in clock 
pulses of l/60th of a second. To get elapsed time in minutes, the clock must be read 
twice and the difference divided by 3600. TIME1 may be replaced by a user's clock read- 
ing subroutine, or it may be removed from the program. 


CONCLUDING REMARKS 

It is not always possible to obtain sufficient detail on some critical parts of the blade 
surfaces by using the TANDEM program. Due to storage limitations on the computer, 
grid spacing may be too large to give the desired detail around small leading- or trailing- 
edge radii or within slot regions. For this reason, a computer program called MAGNFY 
has been written to obtain a solution on a finer mesh in a small part of the blade -to -blade 
region. MAGNFY is described in reference 13. 

After TANDEM was written, it was realized that the TANDEM program was sig- 
nificantly improved over the 2DCP program (ref. 3) for a single unslotted blade. Hence, 
TANDEM was modified to solve the same problem as 2DCF. This modified program, 
called TURBLE, is described in reference 14. The coding in TURBLE is simpler and 
more foolproof than that of 2DCP. Also, TURBLE allows more interior mesh points in 
the solution region, and has its own error package independent of the Lewis computer 
system. It is intended that TURBLE should supersede both the 2DCP and the 2DINCP 
(ref. 11) programs. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 4, 1968, 

126-15-02-31-22. 


97 



APPENDIX A 


FINITE-DIFFERENCE APPROXIMATION 

An approximate numerical solution for the stream function u can be obtained by 
finite-difference methods* These methods involve first establishing a rectangular grid 
of mesh points in the region, as shown in figure 14* Then at each point where the value 
of the stream function is unknown, a finite -difference approximation to equation (1) can 
be written* Adjacent to the boundary, the boundary conditions are included* If there are 
n unknown values, n nonlinear equations are obtained in n unknowns* The equations 
are nonlinear since the coefficients involve the density, which depends on the solution* 

The equations may be solved by an iterative procedure, with two levels of iteration* The 
inner iteration solves a linearized equation, and the outer iteration makes corrections to 
the linearized equation so that the solution converges to the solution of the original non- 
linear equation* 

First, the inlet absolute total density is used for determining the coefficients of the 
finite -difference approximation to equation (1). This results in n linear equations. 

These linear equations may be solved iteratively by successive over relaxation, as de- 
scribed in references 10 and 11* This solution is an approximate solution of equation (1) 
for the stream function* This approximate solution may be differentiated numerically to 
obtain approximate velocities from equations (2) and (3)* The approximate velocities are 
then used to obtain a better approximation to the density at each point, and the coefficients 
of equation (1) are recalculated by using new densities. Thus, the solution to the non- 
linear equation (1) is approached by a sequence of solutions to linear equations* 

A typical mesh point with the numbering used to indicate neighboring mesh points is 
shown in figure 17* The value of the stream function or the other variables at 0 is denoted 


2 


8 h 3 

h 2 

h 4 

3 

0 

h l 

(A8)j - ^L 

1 


» FFl 


Figure 17. - Notation for adjacent mesh points 
and mesh spaces. 
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by using the subscript 0, and similarly for the neighboring points. It can be shown 
(ref. 10) that equation (1) can be approximated by 


r 

[ h l< h l + 


2u 2 2u 0 ' 

+ hg) ,1 2^1 + b 2* h l b 2 




2u a 


2u£ 
4 ) h 3 h 4 


[ h 3^ h 3 + h 4 J + V 

1 /pz^A/^^ y sin g o 

% \ h l + h 2f \ h l + h 2/ . r 0 


(Al) 


Vi’Vsl ( U ± ' u 3\ 
Vo<WJ y 3 + V 


— Vo ata “0 

w 


where h^ = r^AS^ and hg = r^At?^ (since r Q = r j = r 2 ). In setting up equations for 
solution, the coefficients of the u^ in equation (Al) must be calculated. This was done 
by expressing equation (Al) as 


4 

“o “ E Vi + k o 


where 


a 


2 



a 34 : 


h 3 h 4 


b 84 " 


a 0 ' a 12 a 34 

b - P 2~ P 1 

12 Po< h l + h 2> 
b 4 p 4 ~ b 3 P 3 sbl a 0 

Vo (ll 3 + h 4> r 0 






(A2) 
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This equation can be used at all interior mesh points, and for mesh points adjacent to the 
blade surfaces BC, ML, and so forth. 

Along the boundary where the value of u is unknown, the equation will vary. For 
example, along the upstream boundary, 3u/3q is known, and a finite -difference approxi- 
mation to (Su/Sq)^ in equation (4) gives 


u 0 = u 4 



- u 4 + 1*4 


/tan 


(A3) 


Similarly, along the downstream boundary, equation (5) gives 


u 0 ~ u 3 



= 



(A4) 


For the points along AB, equations can be derived by using the periodic boundary con- 
dition. H the point 0 (fig. 18) is on the boundary between A and B, the point 1 is outside 
the boundary. However, it is known that u 1 = u.^ - 1 where the point 1, s is a dis- 

tance s above point 1 in the 6 -direction, as shown in figure 18. Substituting this con- 
dition in equation (A2) gives 


4 

u 0 = a l u l,s + g a i u i- a l + k 0 


<A5) 
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where a^ is the same as defined in equation (A 2). 

The points along MN are not part of the solution regions, since the value of the 
stream function at each of them is just 1 greater than the corresponding point along AB. 
The equation for the first mesh line below NM must be modified, however. In this case 
°2 " u 2 - s + lj where tiie P 0 ^ 2, -s is a distance s below point 2 in the negative 6 - 
direction, as indicated in figure 19. Substituting this condition into equation (A2) gives 

u 0 ~ a l u l + ^“2, -s + a 3 u 3 + a 4 u 4 + a 2 + k 0 

In a similar manner, equations can be derived along the other boundaries (FG, HI, 

CD, and KL,; and DE and JK for the nonoverlapping case, fig. 4) where a periodic condi- 
tion exists. Rather than give the equation for every possible case, it is easier to state 
the rule for modifying equation (A2). E an adjacent mesh point i {fig. 17) is outside the 
mesh region (along a periodic boundary), two changes must be made: 

(1) Change the subscript of u from i to i, s if the periodic boundary is along the 
bottom of the mesh region. Change the i to i, -s along the top of the region. 

(2) Subtract from k^ E the periodic boundary is along the bottom of the mesh 
region. Add to kg along the top of the region. 

One of equations (A2) to (A6) can be applied to each mesh point for which the stream 
function is unknown in the region of interest, giving the same number of equations as there 
are unknowns. These points where the stream function is unknown are referred to simply 
as unknown mesh points. 


N Z M 



<A7) 


101 


This system of n equations is represented in matrix form as 


Au = k 


(A7) 


T 

where u = (u^, . . . , u fl ) is a vector whose components are the unknown values of the 
stream function, A is the coefficient matrix of equations (A2) to (A6), and k = 

(kj, . . . , k n ) T is the vector whose components are the known constants of equa- 
tions (A2) to (A6). If the mesh size is sufficiently small, the coefficients a^ to a^ in 
equation (A2) will all be positive (for any given continuous functions b and p). In this 
case, the coefficient matrix A is irreducibly diagonally dominant, and there is a unique 
solution to equation (A2) {ref. 10). 

The solution to equation (A2) is obtained by using two levels of iteration. The inner 
iteration consists of solving equation (A2) by using fixed values of p based on the pre- 
vious inner iteration. The inner iteration is successive over relaxation using an optimum 
over relaxation factor S2, as described in reference 11 (p. 77). The iterative procedure 
is given by 



for i = 1, 2, . . . , n (A8) 

where £2 is the overrelaxation factor. The a-, are the elements of the matrix A, and 

1J n 

the are the components of the vector k of equation (A7). The u^ are the initial es- 
timates of the and are obtained from the previous inner iteration. 

The outer iteration consists of making corrections to the coefficients so as to finally 
obtain a solution to the nonlinear equation (1). The optimum value of ft can be deter- 
mined as described in reference II (appendix B). The optimum value of ft will vary 
slightly each time the coefficients are corrected; however, the change is usually small, 
and it has been adequate to use the same overrelaxation factor for the entire calculation. 
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APPENDIX B 


NUMERICAL TECHNIQUES USED IN PROGRAM 
Calculation of Velocity and Density 

When the stream function u has been calculated, it is then possible to calculate the 
derivatives 3u/3m and 3u/3f? by numerical techniques. Then, with equations (2) and 
(3), and since W 2 = + W 2 . values for pW can be calculated. It is assumed that the 

values of w, X, r, y, c , T'. , and pV are all fixed and known. Then p, and hence 
pW, is a function of W. The product pW has its maximum value when W = W cr . If 
pW is less than this maximum value, there are two values of W which will give this 
value of pW, one subsonic and the other supersonic. It is desired to find the subsonic 
value of W corresponding to the given value of pW. The method used is Newton's 
method, which converges quadratically. 

It is necessary to express pW as a function of W. The static temperature T may 
be expressed as a function of W and r by {see ref. 15, eq. (3)) 


T[ 


= 1 


W 


2 + 2wX - (o>r) 2 


m 


2c Tl 
pin 


With the assumption of isentropic flow 



and the following equation is obtained: 


(Bl) 


(B2) 


pW =p| n W 


W 2 + 2o)X - («r) 2 


2c Tl 
pm 


1 

7-1 


For Newton's method, the derivative with respect to W is needed, 


(B3) 
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2-y 


d(pW) = 
dW 


W 2 p! 

M in 


yRT| 


in 


x _ W 2 + 2wA - (<ur) 2 

y-i 

+ p ln 

" 

1 _ W 2 + 2cuA - (tor) 2 

2c T! 

2c T! 

pin 


pm 


1 

r-i 


(B4) 


Suppose that (pW) g . y is a given value of pW. A first estimate of W is 


(P w ) giv 

Wn 


(B5) 


m 


Then, using Newton's method. 


W n + l = W n + 


(pw) giv - 

p(W n )W n 

d(pW) 


dW 

w=w„ 


n ® 0, 1, 2, . . . 


(B6) 


Since the convergence is quadratic, only a few iterations are needed, and the relative 
change in W is an excellent measure of the relative error in W R . If an estimate for 
W is available from a previous iteration, this value is used for Wq instead of using 
equation (B5). The algorithm given by equation (B6) is done by subroutine DENSTY. 


Calculation of Prerotation X 


The input information for the program determines the value of A = (rV^)^. The 
average value of (pW)^ can be calculated by 


<pw) Ie = 


w 

r le sb le cos 0 le 


(B7) 


where /3^ g is the average value of f3 across BM. The value of W can be estimated by 
dividing this value of (pW)^ by p!^. Then A can be estimated by 

X = r le (W le sbj ^le + wr le> < B8 ) 


where W^ e is the average value of W across BM. From this a better value of Pj e is 
calculated by 
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1 

nr- 1 


(B9) 


Pie = Pin 


1 - 


W f e + 2uA - (or le )' 

2c T- 

pin 


Use of this value of gives a better estimate of the value of W^ g , and then iteration 
can be used with equations (B8) and (B9) until there is a negligible change in P lg . This 
calculation also gives the value of W^ 0 along BM. These calculations are performed in 
PRECAL. 


Calculation of Critical Relative Velocity W cr 

For reference, the critical relative velocity W 0r is calculated at blade leading and 
trailing edges. This is given by 


W 


2 

cr 


= 2 r R t" 

y + 1 


(BIO) 


where 


T” = T! - 2(t,x ~ ( (Ur )^ 


m 


2c 


(Bll) 


This calculation is performed by PRECAL. 


Calculation of Maximum Value of Mass Flow Parameter pW 

The mass flow parameter pW attains its maximum value when W = W . For ref- 
erence, the maximum values of pW along BM and along FI are computed by the program. 
The maximum value of pW is calculated by 


(pW)_ = p* w „ 
'max in cr 


W cr + 2(J)X “ <" r > 2 
1 - - cr 


1_ 


2c T! 

pm 


(B12) 
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where W is calculated by equations (BIO) and (B 11). 

vi 


Calculatfon of Flow Angles p along AN and GH 

If the radius or stream -channel thickness b is not constant in the meridional direc- 
tion, the free- stream inlet and outlet flow angles 0 change along the meridional axis. 
(By free- stream velocity or flow angle we mean the velocity or angle that would exist at 
a point of the stream channel based on conservation of angular momentum, either up- 
stream or downstream of blade). The following relations hold for free -stream condi- 
tions: 


W, 


tan 0 = 


W 


m 


W 0 - V e - ur 


rVg = Constant 


w m = 


w 


pbrs 


J 


(B13) 


From this we can derive the following equation for the f ree-stream angle 0 at any point 
along the meridional axis, when it is known at some other reference coordinate of 
m = m*. 


tan 0 = 


tan 0 + / p \ “ r 2 )ps 

b* \P*) 


w 


(B14) 


Equation (B14) may be used at either inlet or outlet to calculate 0 in or 0 t . This 
requires iteration, since p is not known until 0 is known. 


Equation for Leading- and Trailing- Edge Radii 

The equation for the leading- and trailing- edge radii is needed. H the radius r were 
constant, 


106 



(m - m*) 2 + r 2 (e - 0*) 2 = R 2 


(B 15) 


where R is the leading- or trailing -edge radius and m*, 0* are the coordinates of the 
center of the radius. Since r changes by a relatively small amount on this circle, it 
was deemed adequate to use this equation with r taken at the leading or trailing edge. 
Equation (B15) is used by the program to calculate coordinates on the leading- and trailing 
edge radii. It is also used to calculate the points of tangency to the spline curves de- 
scribing the rest of the blade surfaces, and to calculate slopes on the leading- and trailing 
edge radii. 


Calculation of Surface Length 


It is often desired to plot the velocities as a function of blade -surface length. For 
convenience, the approximate blade-surface length is calculated by the program. The 
calculation is based on straight-line distances between each vertical grid line on the blade 

surface. Eh. is the spacing between vertical grid lines, r. the radius at the i** 1 verti- 

^ th ■ 

cal grid line, and Q. the coordinate of the i vertical grid line, the surface length S n 

to the n grid line is approximately 


n 




This may be in error near the leading or trailing edge, but is quite accurate over most of 
the blade surface. 
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